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ABSTRACT 


A new numerical method, based on the Vortex Method, for the simulation 
of two-dimensional separated flows, has been developed and tested on a wide 
range of cases. The fluid is incompressible and the Reynolds number is high. 

A rigorous analytical basis for the representation of the Navier-Stokes equa- 
tions in terms of the vorticity is used. It includes an equation for the control 
of circulation around each body, which has sometimes been overlooked. 

In the Vortex Method the vorticity transport equation is solved numerically 
in a Lagrangian reference frame, by following elementary vortices. The 
resulting method is grid-free and concentrates its points in the regions of 
steep gradients; it also allows a simple and exact treatment of the far-field 
conditions. It is well adapted to the modeling of transport phenomena. 

The Vortex Method has been criticized for its handling of the viscous 
effects. In this study most of the effort has been devoted to understanding and 
controlling the parasitic numerical effects, and to reproducing the true physi- 
cal effects. This was achieved by coupling an inviscid outer flow (computed 
by the Vortex Method), with a viscous boundary layer flow (computed by an 
Eulerian method). 

Two significant advantages of this new version of the Vortex Method are 
the capacity to treat bodies of arbitrary shape, and the ability to accurately 
compute the pressure and shear stress at the solid boundary. These two 
quantities reflect the structure of the boundary layer. 

Several versions of the method are presented and applied to various 
problems, most of which had massive separation. The comparison of its 
results with other results, generally experimental, demonstrates the reliability 
and the general accuracy of the new method, with little dependence on em- 
pirical parameters. Many of the complex features of the flow past a circular 
cylinder, over a wide range of Reynolds numbers, are correctly reproduced. 

The method appears to incorporate many of the physical mechanisms of 
separated flows, and the dependence on Reynolds number has been obtained. 
Its accuracy, when experimental results are taken as a reference, is limited 
mostly by difficulties in modeling turbulence, and by the two-dimensional 
assumption. 
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Nomenclature. 


Roman symbols. 


(a, b) A set of Lagrangian coordinates. 

A m Area of solid region S m 
Bf Intensity of the buffer vortex sheet, 
c Chord of airfou. 

Cd Drag coefficient: Cd — dras'/(.5/)c 2 f7 0O 2 ) 

Ci Lift coefficient. 

C n Normal force coefficient. 

C m Moment coefficient: C m — moment / (.5 pc 3 f/oo 2 ) 

C p Pressure coefficient: C p = (p — Poo)/{-5pUoo 2 ) 

Cpo> Cpi«» (7pi c Fourier coefficients of the pressure during the dynamic stall. 
D 0 Parameter in merging device. 

Di Distance from i th vortex to the wall. 
d /dt Eulerian time derivative (fixed point in space). 

D /Dt Lagrangian time derivative (particle). 

F Fluid region. 

i Imaginary complex number: t 2 = — 1 
i,j Indices of two vortices. 

k Reduced frequency of the dynamic stall pitching. 

I Side of the square cells used for the taylor expansions. 

K,L Indices of two cells for Taylor expansions. 

L Length scale associated with the solid body, 

m Index of a solid. 

M Number of solids. 

n Coordinate normal to the wall in inner region, 
n Unit vector normal to dS. 
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n Number of cells for approximating the vortex interactions. 
N g Number of grid points in a Finite Difference simulation. 

N v Number of vortices. 

N w Number of points along the wall. 

p Pressure (divided by the density), 

r Vector representation of (x,y): 



R (large) distance from the origin in an asymptotic expansion. 
Re Reynolds number. Re = \Uoo\L/i/ 

R 0 Distance from the creation points to the wall. 
s Coordinate parallel to the wall in inner region. 

As Size of intervals along the wall. 

S Solid region (union of the S^’s) 
dS Boundary of S (union of the dS^'s). 

S m Interior of the m th solid. 
dS m Boundary of S m 

t Unit vector tangent to the boundary dS 

t Time. 

u, v Velocity in x and y directions. 

U Velocity vector: 

u =0 

U m (z>y) Velocity of the solid material in the m th body at ( x,y ). 
U m o Reference velocity of the m th body (at z = y = 0) 

Uoo Velocity at large distances (in general aligned with x axis). 
Vo Tolerance in the merging device. 
x, y Coordinates. 

z Complex representation of ( x, y) : z — (x -f iy) 

Z Complex variable, used in model equations. 
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Greek symbols. 


f incidence of the airfoil. 

<* 0,07 mean and amplitude of a during the dynamic stall. 

0 intermittency factor for the turbulence model. 

7 normalized core vorticity distribution. 

T Circulation: closed line integral of the velocity. 

1\ Circulation of ith vortex. 

8 thickness of the numerical viscous region. 

8 centroid of the inner region vorticity in the n direction. 

6 filtered version of 5*. 

A Laplace’s operator: A = d x * + <9 y a 

At Time step of the numerical integration. 
e Amount of artificial dissipation. 

rj Function of regularization of the velocity, associated with 7. 
V Gradient operator: 



H Coefficient of viscosity. 
u Coefficient of kinematic viscosity: u = n/p. 
p Density of the fluid (p is constant and omitted in general) 

u Vorticity. 

fi m Angular velocity of m th solid. 

¥ Stream function, 
u Core radius. 

r Viscous shear stress at the solid wall. 


Subscripts and superscripts. 


00 value at large distances (freestream). 
z complex conjugate of z. 
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I) INTRODUCTION. 


1) Description of the proL>*^. 

Separated flows, in which the fluid fails to smoothly follow the solid sur- 
face, in contrast with attached flows, are generally more complex and more 
challenging to measure or predict. It was shown by Prandtl that, for the 
same conditions, the viscous equations and the inviscid equations can have 
very different solutions even if viscosity is extremely weak, precisely because 
the viscous flow might separate while the inviscid flow does not [1] . Flows 
at high Reynolds numbers (low viscosity) are thus very sensitive and, by and 
large, we only have a qualitative and sometimes simplistic knowledge of their 
behavior. 

In most designs, separation is undesirable since it results in inefficient 
operation, with high drag or loss of pressure, or even it leads to a dangerous 
situation like stall. However many devices, like wings or diffusers, often 
operate on the verge of separation. In other cases separation is present in the 
design conditions, where it is undesirable but unavoidable, like on aircraft 
tails or cars, or a normal feature of the flow, like on a three dimensional 
wing. 

A relevant example is the flow around the retreating blade of a helicopter 
in translation [2] , (3] . The blade experiences large and rapid changes of 
incidence and velocity, sufficient to cause stall with strong unsteady effects. 
The blade may even move with the trailing edge forward for part of the 
cycle, which always causes separation. The dynamic stability of the system 
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strongly depends on the aerodynamic loads, especially the pitching moment, 
and these exhibit very significant non-linear and hysteretic effects. Figure 
1 illustrates bow much the loads can differ during dynamic stall from static 
loads at the same incidence (Fig. 1 is a personal communication by W. J. 
McCroskey. The data appeared in [60]). Accurately predicting dynamic stall 
woul<' thus be very useful, and should be possible in the near future (at least 
for two-dimensional flow), thanks to improvements in numerical methods and 
computers. 

Experimental results demonstrate the strong sensitivity of separated flows 
to details of body geometry, for instance surface roughness [4] , and of course 
to the Reynolds number [1] . The best known example is the circular cylinder: 
this flow still exhibits dramatic changes at Reynolds numbers of one million 
(Fig. 2 and 3). Wind tunnel tests are less reliable when fine viscous effects 
are involved than when the compressibility effects dominate, because the 
Reynolds number depends on model size, while the Mach Dumber does noi, 
s t d because separation is influenced by wind-tunnel turbulence. Analysis 
alone has not been able to produce many results, mainly because separated 
flows can rarely be treated as slightly perturbed from a known exact solution: 
they are not very accessible to small disturbance theories. Free streamline 
theory made use of the observation that, in many cases, drag depends mostly 
on the forebody shape (upstream of the separation point) and very little on 
the part of the body which is inside the wake [5] . Pressure also appears to 
be almost constant in that . ^gion. The idea developed was to treat the wake 
as a B dead-watei" region and to assume a constant pressure in that region; in 
general the value of the base pressure is determined empirically. This theory 
did produce some good results [6] , but a method that ignores the unsteady 
character of the flow cannot be expected to be very accurate; it relies very 
much on empiricism. Thus there is an existing need for the development of 
numerical methods capable of solving either the full Navier-Stokes equations 
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or a high level approximation to them without relying on much empirical 
input. Once these methods have reached an acceptable level of accuracy, 
they are expected to be much faster and less expensive than large scale wind 
tunnel tests, and may be more accurate if they are carefully validated against 
flight tests. 

Whereas accurate and practical numerical methods are available to com- 
pute attached flows [7] , similar methods do not exist for separated flows, 
which are vortical. The slowly-varying, attached, irrotational flows are very 
amenable to finite difference or finite element methods, and to an Eulerian 
formulation. Some of these methods can also treat separated flows, but 
obtaining accurate results becomes extremely costly at moderate or high 
Reynolds numbers [8] . One alternative is the Lagrangian "Vortex Method” 
[9] . This method provides a description which is better adapted to high- 
Reynolds number, vorticity-dominated, unsteady flows, and should result in 
a greater accuracy for a given level of computing resources. 

The first objective of this work was to study the capabilities of the Vortex 
Method, review its inherent strengths and weaknesses, especially in the con- 
text of two-dimensional separated flows, and remedy some of the weaknesses. 
The other objective was to develop a reliable and accurate computer program, 
based on the Vortex Method, for the simulation of a general class of separated 
flows. This program has been validated by systematic comparison with known 
results, and is beginning to be used as an active research tool to investigate 
some candidate designs, in parallel with wind tunnel tests. 

The flows to be considered are viscous flows past two-dimeasional solid 
bodies in a uniform stream. Only incompressible flows are considered. I he 
incompressibility limitation is associated with the Vortex Method. The two- 
dimensional restriction is not, but simai. ,ing two dimensional flows is a first 
step ani reflects the "state of the art”. (The extension to three dimensions 
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would not be straightforward, but it is certainly possible [9].) We consider 
here one or more bodies aM they may be in non-uniform motion. Even if the 
motion is uniform, the flow is likely to be unsteady with a possible periodic 
character. Frequently separation of the boundary layer will occur as a result of 
the body being bluff or at high angles of attack. Large vortical structures will 
appear and form a wake having a turbulent character, md these structures 
will strongly influence the load’s on the body. Their subsequent decay in the 
wake far downstream is of less interest because of their small influence on the 
loads. Again, typical examples are the flow past a circular cylinder, and the 
static or dynamic stall of au airfoil. 

2) Equations. 

The Dv havior of isotropic viscous fluids is described by the Navier-Stokes 
system of partial differential equations. The independent variables are the 
cartesian coordinates x and y and the time t. In the conventional formulation 
the dependent variables are the velocity vector U = (u, v) and the pressure 
p. The density p is constant since the fluid is incompressible; the symbol p 
will actually be taken to represent the ratio p/ p, and p will be omitted in 
the writing. Similarly, the coefficient of viscosity p is divided by p to yield 
the kinematic viscosity v. The dependent variables are defined in the fluid 
region, that is the region of the plane exterior to the solid. Since the fluid is 
incompressible the problem involves only U and p, and the following system 
of equations prevails: 

(continuity) V.U = 0 (1) 

QTT 

(momeutum) — — f- U.VU = — Vp -f- i/AU (2) 

at 

where V is the gradient operator and A is Laplace’s operator. 
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The boundary conditions are as follows. At large distances from the body 
U tends to the "freestream velocity” I/,*,. At the boundary with a solid the 
velocity U of the fluid equals the velocity of the solid material. No boundary 
conditions are needed for the pressure. Initial conditions for U at time 0 are 
also considered given. 

Instead of the conventional ” (u, v, p)” formulation, the vorticity formulation 
is sometimes useful. The role of vorticity in the dynamics of the problem 
considered here is crucial, and a more efficient method is likely to result if 
the vorticity is treated directly. It is defined by: 


_ dv du 

W dx dy 


(3) 


In two dimensions a; is a scalar quantity and is interpreted as the local angular 
velocity of the fluid (multiplied by 2). 

We will show that, owing to the boundary conditions imposed on U, there 
is a one-to-one correspondence between an incompressible velocity field U and 
a vorticity field u. This allows one to develop a solution by focusing on the 
vorticity. 

The vorticity obeys a well-known conservation law. Taking the curl of 
Equation 2 and using Equation 1 we obtain: 


Doj 

~Dt 


— — I- U.Vw = u&u) 
dt 


(4) 


Equation 4 describes how vorticity is convected by the velocity field and 
diffused by viscosity. In two dimensions there is no term corresponding to 
"vortex stretching". Thus Equation 4 is of the same type as the equation 
governing dye concentration. If dye is released by the solid it stays in 
streaks that trail the solid and are confined to the wake. So will voriicity. 
The difference between dye and vorticity is that dye concentration does not 
interact with the velocity (it is a "passive" scalar) while vorticity and velocity 
are related by Eq. 3. 
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The incompressibility condition is now implicit and the pressure term drops 
from the equations. And as a result the number of unknowns is reduced 
from three to one. The main difficulty is in deriving appropriate boundary 
conditions (or conditions of another type) to form a complete system with 
Eq. 4. Many approaches exist in the literature and the one taken for this 
study will be described in detail later. 

We can now turn our attention to the principal numerical methods that 
have been proposed to solve either Eqs. (1,2) or Eqs. (1,3,4). 

3) Related investigations. 

The finite difference method is the prevailing method in Computational 
Fluid Dynamics, as opposed to finite element methods, spectral methods, 
vortex methods, etc, and will be reviewed first. The finite element methods 
will not be described. They are quite similar to the finite difference methods 
and are receiving more and more attention because they are f vmally more 
accurate. They ai-e probably less mature and certainly less widespread, at 
least in the English literature on fluid mechanics. Spectral methods can be 
extremely accurate, but are still much less versatile than the other methods; 
they have been used only with very simple geometries (periodic flows, or 
channels) and not for flows around solids. For that reason, they too will not 
be described. Finally, the Vortex Method is a promisiDj although not very 
mature alternative to finite differences for the simulation of incompressible 
vortical flows. The method will be introduced and its literature reviewed 
after the finite difference methods have been considered. Then the relative 
advantages of the two methods will be assessed. 

a) Finite difference methods. 

The finite difference method is very well known [10] . Out of the large 
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number of finite difference publications, we shall describe only a few out- 
standing studies, whhh seem to be capable of treating two-dimensional flows, 
with large separation, at Reynolds numbers of at least 10 3 . As a rule, the 
difficulty increases with the Reynolds number. 

Ail of these studies used the Eulerian frame of reference and sotoed either 
Eqs. (1,2) or Eqs. (1,3,4) by finite difference approximations on a grid 
that does not evolve in time. With either formulation there is a variety of 
finite difference schemes available, time advance schemes, boundary condition 
procedures, and turbulence models if applicable. 

In 1961 Thoman and Szewczyk treated the flow over a circular cylinder 
for Reynolds numbers ranging from 1 to 3. X 10 5 [11] . They used two 
overlapping grids: one near the surface and an outer grid, extending to only 
5 diameters. Freestream conditions were imposed over most of the outer 
boundary. It is difficult to estimate a priori how much this affects the solution, 
compared to a situation in which the disturbances are allowed to extend much 
farther than 5 diameters. Thoman and Szewczyk used an upwind differeLce 
scheme to stabilize the computation. They recognize the important fact 
that this scheme is dissipative, in the sense that its stability comes from 
a numerical dissipation of the energy, not the physical dissipation. (The 
elementary form of upwind differencing introduces enough diffusion to bring 
the effective cell Reynolds number down to 2.) Thoman and Szewczyk carried 
their computations up to the onset of the drag crisis and the average drag 
they found was very accurate. They did not report results at higher Reynolds 
numbers. The pressure distribution was accurate up to a Reynolds number 
of 400 and quite inaccurate at 3. x 10 5 , although fortuitously the drag did 
not reflect it. The Strouhal number was too low by about 30%. The results 
were quite good, but the accuracy of the upwind scheme was a subject of 
controversy. 

Ten years later Jordan and Fromm treated the circular cylinder for 
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Reynolds numbers ranging from 100 to 1000 [12] . They used a log-polar 
grid of large extent (187 diameters) and an outer edge condition devised to 
allow the solution to oscillate freely. The time history of lift, drag, and 
torque clearly showed that a limit cycle was reached. The drag and shedding 
frequency were accurate but again the pressure distribution was not as satis- 
factory. The authors estimated that their computations should be considered 
as accurate up to Re — 400. 

In 1977 Mehta computed the dynamic stall of an airfoil [8] . He used the 
vorticity formulation, Eqs. (1,3,4), a conformal mapping from the airfoil 
to a circle, and finite difference approximations. The numerical boundary 
conditions imposed at the outer edge were chosen to constraint the solution as 
little as possible. A very elaborate implicit program was used, to obtain high 
order accuracy and reasonable running times. Implicit time marching schemes 
are more stable, numerically, than explicit ones. The flow was incompressible 
and the simulation "direct” (no turbulence model) with Reynolds numbers 
up to 10 4 considered. Good agreement with flow visualizations was obtained. 
All the qualitative features of the flow were reproduced, but quantitative 
comparisons were not reported. 

Wu treated the flow around a circular cylinder and around a stalled air- 
foil by an original method [13] . Wu presented a very good description and 
justification of his vorticity formulation in [14] . For the numerical method, 
he computed the vorticity on a grid, but only the cells that contained vor- 
ticity were active. This helped reduce the number of points, like in the Vortex 
Method, except that here an active cell could never be passive again: the 
computational domain could only grow with time. Also, while the vortical 
domain is formally infinite (because of viscosity), Wu kept it finite by activat- 
ing a cell only if it contained more than an arbitrary "low” level of vorticity. 
The irregular domain is expected to make the vectorization of the program 
difficult. The velocity at the grid points was computed by Biot-Savart in- 


tegration and the vorticity equation was solved in Eulerian coordinates by an 
explicit method. In a recent paper Wu and Gulcat treat separately the wake, 
the irrotational region and the attached boundary layer [15] . By adopting 
the simplest possible level of description for each region, they save significant 
computer time. Wu obtained very accurate results for the cylinder at low 
Reynolds numbers. At a Reynolds number of 4. x 10 4 Wu and Gulcat ob- 
tain what appears to be a very good pressure distribution and a good drag 
coefficient. However they compare the experimental pressure, averaged over 
a long time, with the computed instantaneous pressure at time 4.8 (based 
on velocity and radius). After such a short time the flow surely has not 
reached its limit cycle. Thus the agreement might be fortuitous. In general, 
Wu produced some very good ideas tut did not always support them with 
sufficient numerical evidence. 

In 1981 Tassa and Sankar treated dynamic stall in compressible turbulent 
flow [16] . They used an implicit finite difference program and an algebraic 
turbulence model. The overall quality of their results was comparable to 
the quality of the results to be reported in the present work. The agree- 
ment between different simulations or different experiments, for this difficult 
problem, is only qualitative. Shocks were not mentioned although the Mach 
number was 0.6 and high incidence angles were reached. It also seems that 
the downstream boundary condition used would not allow circulation to leave 
the computational domain; this is a problem with any method that solves the 
equations on a finite domain. 

The study by Shang, in 1982, treated compressible flow around a circular 
cylinder [17] . Shang plans to extend it to three dimensions. Accordingly, 
he used the primitive variables (density, velocity, pressure, energy) and the 
compressible equivalent of the system of equations (1,2). The computational 
domain extended to 30 diameters; ” non reflecting” boundary conditions were 
applied at the outer edge to minimize the constraint introduced by the finite 
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domain. The explicit McCormack scheme was used. The program was fully 
vectorized on the CRAY computer. The Mach number was 0.6 and the 
Reynolds number 1.67 x 10 6 , which is quite high, but no turbulence model 
was implemented. This suggests that the computation was stabilized by 
the numerical dissipation of the McCormack scheme, which might be much 
stronger than the molecular dissipation, depending on the grid. The drag and 
shedding frequency agreed very well with experiments. The average pressures 
were not reported. The lift exhibited a markedly non-harmonic behavior; this 
is not mentioned in textbooks, but is consistently observed in experiments 
and in computations, both finite difference and vortex. 

In 1982 Davis and Moore treated the incompressible flow past rectangles at 
Reynolds numbers between 100 and 2800 [18] . In that range, the molecular 
viscosity still has a significant effect and the flow characteristics depend on 
the Reynolds number. The finite difference scheme was chosen to provide a 
smooth solution with a minimum of numerical dissipation. The freestream 
conditions were imposed at the outer boundary, except on the downstream 
face where the numerical boundary condition was chosen to allow vortices to 
cross the boundary. The grid was adapted to the rectangular shape and it 
might be difficult to extend the program to other shapes. Satisfying agree- 
ment with experiments was obtained, especially at low Reynolds numbers. 
They estimated that the upper limit for good accuracy was about 1000. 
Computations at a higher Reynolds number will require a very fine grid, or 
a turbulence model, or both. A remarkable feature was that, while the flow 
at Re — 250 was very regular, with the lift signal almost a pure sinusoidal 
function of time, at Re = 1000 the shedding was much more irregular. 

b) Yortex methods. 

Before a discussion of the literature is presented, the basic idea of the 
Vortex Method will be introduced and its most salient features discussed. 
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The Vortex Method was designed in an attempt to provide a more natural 
and efficient description of the eddies and of the vorticity they carry. The 
method represents the vorticity field as the sum of a large number N v of 
mobile functions with small supports: 

= r< ^l r_r »D ( 5 ) 

*=i 


where r* = (z,-, jk) is the center and I\ is the circulation of the i th vortex, 
and 7 is the function of regularization or ” core shape” . 7 is smooth, has a 
small support, and an integral of 1. ~’ypically, 7 is a Gaussian. This provides 
a very convenient description of the vorticity. The main advantage, when 
external flows are treated, is that in the large irrotational region no vortices 
will be needed. This saves large amounts of memory and allows vortices to 
be concentrated in the wake, where resolution is needed. 

Dynamically, these blobs follow the fluid, like particles. This is a 
Lagrangian description. They retain their circulation in time, so that total 
vorticity is conserved; this corresponds to the inviscid fluid equations 


dTi 

dt 


= 0 


( 6 ) 


|f = W(K. <) (T) 

Equations (5,6,7) give a closed problem involving only the r^’s and IYs, 
provided that U can be calculated from u. U needs to be known only at 
the vortex locations and not in the irrotational region; with incompressible 
fluid, this can be achieved by application of the Biot-Savart lav. On the other 
hand, the main disadvantage with using the Biot-Savart law is that it makes 
each vortex interact with all the other vortices at every time step, which is 
very costly. 
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The advantage of a Lagrangian description for the solution of the inviscid 
version of Eq. 4 is obvious: in Lagrangian coordinates u is constant in time. 
The transport of any quantity is always treated better by having the quantity 
travel across the domain rather than by transferring the quantity from a fixed 
grid point to the neighboring points. As a result, the Vortex Method has no 
obvious numerical dispersion and possibly less numerical diffusion than an 
Eulerian method (this last point will be discussed in more detail). The Vortex 
Method also turns out to be much more stable than most Eulerian methods. 
Large time steps can be taken as long as the accuracy is sufficient. 

The efficiency of the Vortex Method, compared to a ”u, v,p ” formulation, 
arises in particular from the exploitation of two assumptions: the fluid is 
incompressible and inviscid. 

The incompressibility restriction is clearly necessary to the Vortex Method 
in its present form (the Biot-Savart law depends on it). With air, it means 
that only flows at low Mach numbers can be treated, such as the flow around 
a landing airplane, or around the retreating blade of a helicopter. Even then, 
high subsonic Mach numbers can appear locally for freestream Mach numbers 
as low as 0.2. So far these effects have had to be neglected. For flows of 
liquids, incompressibility is obviously a good assumption. 

The inviscid restriction is more controversial. The convergence of the 
algorithm to the solution of the inviscid equation has been mathematically 
demonstrated (in the absence of boundaries) [19] . Explicitely adding the 
viscous term v&u is not convenient in a Lagrangian reference frame because 
it involves derivatives with respect to the Eulerian coordinates. On the other 
hand the method often reproduces viscous behavior, especially around solids, 
even though it is based on the inviscid equation. For years this feature has 
been used to simulate viscous flows with an "inviscid” method. This will 
be made clearer by use of some theoretical arguments and some numerical 
experiments described in this report. 



The method has been studied and refined for decades, without becoming 
operational and widely used. This is partly true because of the viscosity issue 
and partly because the problems it is applied to are very difficult for any 
method. We shall limit our review to papers treating flows past solids. Good 
research has been done on flows without boundaries, but the difficulties these 
simulations raise are quite different: in these cases viscosity actually plays a 
negligible role. 

Bryson, in 1959, used a very simple model to represent the flow around 
a circular cylinder, with one pair of symmetrically placed vortices which 
moved away from the cylinder and gathered circulation with time [20] (Fig. 
4a) (Bryson did not use Eqs. 6 and 7). Thus, flow separation was assumed 
but viscosity was not accounted for otherwise. This rather empirical model 
served well for a short time after a rapid start. It was intended for use in 
a slender body analogy: the steady three dimensional flow past a slender 
cone at angle of attack is analogous, cross-section by cross-section, to the 
two dimensional time-developing flow past an expanding circle in translation. 
The two flows have many common features, including the formation of two 
symmetric vortices, followed by a loss of stability and an asymmetric state 
with a side force. This side force can affect the control of airplanes with long 
noses. 

At the next level of complexity, a large number of vortices are used and 
follow the fluid (Eqs. 6 and 7 are applied) and symmetry is not imposed 
[21] , [22] , [23] . New vortices are added at the separation points, which 
are either obvious (a corner on the body) or known empirically (the leading 
edge of an airfoil or an assumed separation point on its top surface, the 84° 
point cn a circular cylinder in the subcritical range, etc ). The strength and 
exact position of the new vortices are chosen in accordance with Prandtl’s 
rule and a so-called "Kutta condition”. Prandtl’s rule states that the flux of 
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separating vorticity is u 2 /2, where u is the velocity at the outer edge of the 
boundary layer. The "Kutta condition” is applied even though the wall is 
smooth; it states that the velocity at the wall, under the separating boundary 
layer, is zero. It would be described better as a selective application of the 
no-slip condition. Thus the boundary conditions are neither truly inviscid 
nor truly viscous. The distinction that is made between " under" and ” over” 
the boundary layer is not very clear, especially when the upstream part of 
this boundary layer is not represented. This is a major source of uncertainty; 
the exact points where the velocity is sampled are quite arbitrary and have 
a strong influence on the results [24] . A mapping from the body shape to 
a circle is used, in conjunction with image vortices, so that the tangency 
condition is satisfied. This is the traditional way to treat inviscid flows. It is 
not as weil adapted to viscous flows, and the method to be described in this 
report actually does not use images. In addition to degrading the accuracy 
(the interaction of a vortex with its image becomes very inaccurate when they 
are close to the w U), the use of mappings and images ' extremely unwieldy: 
accurate conformal mappings for arbitrary shapes are not readily available. 
The viscous "no-slip” condition is satisfied only where the Kutta condition is 
applied. 

Some authors allow the vortices to emanate from a fixed point on an 
ordered shear layer (Fig. 4b). In general this requires a redistribution of the 
vortices at each step to keep the curve smooth [25] . From time to time, the 
shear layer is cut on an empirical basis to allow the formation of the Karman 
street [22] . Other authors do not link the vortices and let them become a 
"jungle” (Fig. 4c). 

In most of the papers of this period, the flow is called "inviscid” and the 
question of how any vorticity can leave the wall is not addressed. The value 
of jeflkient of viscosity is irrelevant. The method tends to overpredict 
the drag and an empirical suppression of vorticity is often used to decrease 



the drag to the desired value [22] ,[26] , [27] . This suppression of vorticity, 
which violates the two-dimensional vorticity conservation law, is sometimes 
presented as a way to account for three-dimensional effects. 

Defifenbaugh and Marshall attempted to couple the Vortex Method to a 
boundary layer, using an integral method for the boundary layer [28] . They 
encountered difficulties in locating the separation point , possibly because of 
the inaccuracies associated with -he release of a single vortex at separation, 
and because of the questionable validity of Bernoulli’s equation in a vortical 
flow. They used a merging device for adjacent vortices. They treated the 
circular cylinder at subcritical Reynolds numbers and concluded that the 
coupling algorithm still had to be refined. They also found that they had 
to destroy some of the vorticity, otherwise the drag came out too large. 
Deffenbaugh and Shivananda proposed a method to treat compressible flow 
at low Mach numbers [26] . Apparently their first attempt was not carried 
further. 

A more ambitious approach was taken by Chorin [29] . A random walk 
displacement is added to the motion of the vortices. This random walk 
reproduces the effects of viscous diffusion statistically. This algorithm treats 
the whole flow as viscous, the Reynolds number is well defined and finite. 
The no-slip boundary condition is used, vortices are present all along the 
wall and the separation of the boundary layer is spontaneous (Fig. 4d). No 
empirical input is needed and the method can now solve problems on its 
own, provided that the resolution is fine enough for the statistical argument 
to hold. Unfortunately, it seems that this would require a huge number of 
vortices and an extremely accurate integration of the transport term, so that 
the scattering it creates does not dominate the random walk [30] , [31] ; the 
random walk idea is attractive but not so practical. 

Chorin treated the boundary condition at the wall in two stages, using both 
sources and vortices. He applied the boundary conditions in collocation form, 
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which is not optimal. The integral form is more costly, since it requires the 
evaluation of logarithms (or arctangents), but it is much more accurate and 
iess sensitive to the non-physical parameters, for example the core radius. 
Chorin treated the circular cylinder; his results for drag are accurate at Re 
= 100, but then seem to decrease monotonically above 100, which is not 
correct physically. With use of the sources Chorin did not need to employ a 
conformal mapping. However both he and his student, Cheer, later returned 
to the image-and-mapping method (probably so that the normal velocity 
would be identically zero at the wall, instead of oscillating near zero). 

Chorin subsequently introduced the "Vortex Sheet Method” in order to 
take into account the widely different scales in the s and n directions c” 
the boundary layer and to reduce the scattering in the direction normal to 
the wall [32] . The region exterior to the boundary layer is treated by the 
"isotropic” Vortex Blob Method with an exchange of vortex elements, sheets 
becoming blobs and vice versa. The circular cylinder and a Joukovsky airfoil 
were treated b 1 , Cheer with this hybrid method [33] . She reports good values 
for the drag of the cylinder (at subcritical Reynolds numbers), but the results 
were not very detailed and the runs seemed to be very short. Chorin and 
Cheer did not use a merging device and had to stop their simulations after 
a fairly short time to keep computer cost under control. In this case, like in 
Wu and Gulcat’s case, such short simulations are questionable as they clearly 
do not reach a true asymptotic state. Although his work left room for many 
improvements it is clear that Chorin showed the way towards a method which 
is powerful and mathematical in spirit, rather than empirical. 

More recently, Lewis independently introduced an image-free form of the 
Vortex Method which is similar to the one that will be presented here [34] . 
It is not clear whether Lewis correctly applied Eq. 10.5 (see below) or an 
equivalent condition. Lewis made use of the advantage of not needing a 
mapping to treat various shapes. Using a modest computer, Lewis introduced 



only a small number of vortices, at only two separation points. The Doundary 
layer war not treated separately. He applied Prandtl’s ruie, although he 
recognized that it is very delicate to apply. Porthouse aed Lewis subsequently 
published results using a random walk model to account for viscosity [35]; 
these results seem to confirm that very many vortices and very short step- 
would be needed for the random walk effect to be meaningful at practical 
values of the Reynolds number. 

c) Finite difference versus Vortex methods. 

It appears that presently neither method ;one can treat high Reynolds 
number flows with the level of accuracy that is needed for engineering. For 
example, reliable quantitative predictions have not been obtained for dynamic 
stall and computing these flows is a matter of research, not of production. 
In some cases, only qualitative agreement is obtained, for instance agreement 
with visualizations. In other cases, the quantitative agreement is good but 
limited to a few numbers l ie drag or shedding frequency. One reason 
is that quantitative and veriLed experimental data are often not available 
for separated flows; these seem to be as hard to measure as they are to 
compute. In addition, the numerical method is still two-dimensional and 
the experiments, even when the geometry is two-dimensional, often ha/e 
significant three dimensional effects. The comparison with exact solutions 
would be a more rigorous test of the accuracy of the numerical results; 
unfortunately, almost no exact solutions are known for separated flows. 

If the Reynolds number is quite low, less than 10 3 , the finite difference 
methods work well, because the solution is very smooth At the Reynolds 
numbers of aeronautical interest, which are of 10® or more, the vortical 
structures in the wake become so small that a very fine grid is needed, 
which requires a very large memory and very sh-'rt time steps. If the grid 
is too coarse numerril diffusion and dispersion can easily dominate physical 
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diffusion. If this happens the simulation is not truly viscous; such a situation 
is often considered acceptable far from the walls, but not close to the wall. 
Reaching higher Reynolds numbers will mostly be a matter of computer 
power, both in terms of memory and speed, and of turbulence modeling. 

The Vortex methods suffer from their imperfect viscous modeling and a 
certain lack of credibility, independent of the Reynolds number. Nearly all 
the studies treated the same shapes: cylinders, ellipses, Joukovsky airfoils. 
Also, too much empirical input was needed. On the other band, the vortex 
programs needed relatively little memory and some of them ran very fast. 
Improving them is more a matter of improving the algorithm. 

A basic advantage of the finite difference methods is that they rest on a 
well established theory of stability and convergence (at least for bounded 
or periodic geometries; infinite domains are not treated in a fully satisfying 
manner). The same cannot be said of the Vortex methods when viscosity and 
boundaries are involved. 

Another advantage of finite difference methods over Vortex methods is that 
they can be extended to compressible flows without major changes. Treating 
compressible flows, with a Mach number above about 0.1, is even easier in 
some cases because it makes the celerity of the signals smaller. The extension 
to three dimensions is also simpler, conceptually, than for the Vortex Method. 

The most significant difference is that the finite difference methods include 
the viscous terms, while the Vortex Method is essentially 5 ^viscid. However, 
the finite difference grid often is too coarse to resolve these viscous terms 
except close to the wall [35] . This effectively removes the laminar viscosity, 
and the energy is controlled by some form of numerical dissipation instead. 
The true advantages of the finite difference method, even over a Vortex 
Method coupled to a boundary layer, are that boundary layer assumptions are 
not involved (therefore no singularities are expected) and that the transition 



from a viscous treatment (fine grid, near the wall) to an effectively inviscid 
treatment (coarse grid, away from the wall) is smooth. 

Another issue is the modeling of the turbulent stresses in the wake. 
Whereas modeling these stresses is reasonably easy in the boundary layer, 
modeling them in the wake is extremely difficult. Thus, some finite difference 
methods include turbulent stresses, but these are evaluated with so much un- 
certainty that the benefit is not obvious in terms of accuracy. The turbulent 
stresses conveniently improve the stability of the finite difference computa- 
tions; this is not an issue with the Vortex Method which is very stable. 

The Vortex Method has only treated bodies in a uniform freestream flow. 
It is planned to extend it to bodies in a uniform shear flow, which will be 
quite simple. On the other hand it would be much more difficult to treat non- 
uniform incoming shear flow. In that domain the finite difference methods 
are still more versatile. 

The main advantages of the Vortex Method are its accuracy in treating the 
convection terms, and th 1 absence of a grid. Generating grids around com- 
plex shapes is not easy [36j , and unless the grid is very smooth the accuracy 
suffers. Furthermore, for many finite difference programs the computational 
efficiency depends on mapping the physical domain to a rectangular com- 
putational domain. Thus, treating several bodies either involves the use of a 
highly distorted grid, or a rather delicate zonal approach [37] ,[38] . These 
difficulties are totally absent in the Vortex Method, at least in its most recent 
versions, which makes it especially attractive for multiple bodies. 

Another advantage is that the Vortex Method effectively includes the 
infinite domain whereas the finite difference methods include only a finite 
domain and require artificial boundary conditions at a finite distance from the 
body. Choosing these conditions is delicate: there is a danger of constraining 
the solution in a hidden way. The Vortex methods need less empiricism in 
this regard. It is also possible to add wind-tunnel wall effects to the Vortex 
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Method (this is quite simple if the flow does not separate from the tunnel 
walls). 

The Vortex Method previously lacked versatility: with the use of conformal 
mappings it was awkward to treat shapes other than ellipses or Joukovsky 
airfoils. The situation has now reversed itself since, as we shall see, recent 
versions of the Vortex Method easily treat arbitrary shapes while avoiding 
grid generation [34] , [39] ,[40] . 

It is not easy to assess the relative computer costs for the two methods. 
In the Vortex Method, interactions have to fcj computed at each time 
step, where N v is the number of vortices. In contrast, many finite difference 
methods require only of the order of N g operations, where N g is the number 
of grid points. This is true for most explicit methods and for the implicit 
methods that have a suitable ordering of the grid. Since both methods, 
in their widely used forms, are second order accurate, the finite difference 
method seems to have the advantage. However in practical situations N v and 
N g are limited and the relevant question is: \ iiich values of N v and N g would 
achieve the desired level of accuracy? Then the memory requirements and 
the running time', could be compared. 

Only experience can answer the question, but two general rules apply. 
First, the Vortex Method will be more competitive if the vortical region is 
small, which makes N v much smaller than N g . Typically, the Vortex Method 
works well with an external flow, but not as well with an internal flow which 
might be filled with vorticity, and thus make N v and N g about equal. Second, 
the Vortex Method nearly always requires less memory, while the running 
times can differ greatly in one sense or the other. Many researchers reported 
extremely short running times for Vortex computation^, but their resolution 
was very coarse and their accuracy questionable. Other Vortex computations 
required hundreds of hours of computer time. 

To allow an evaluation of the method used in this study, the run time used 
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by the computations will be reported in the "Results” section. 

4) Summary of the evolution of the present method. 

The starting point of this study was an algorithm written by R. Rogallo 
(NASA Ames C. F. D. Branch, unpublished work). It was similar to Chorin’s 
1973 method [29] , but it included a merging device and made use of the 
integral form of the boundary condition, which is more efficient. It also had 
images. During the present study, which was also done at NASA Ames, 
the algorithm underwent three mutations, resulting in the versions KPD 1, 
KPD2 and KPD3. 

KPDl makes use of the new boundary condition (without sources, images 
or conformal mapping), but does not employ a boundary layer. It is versatile, 
robust, and accurate for flows that are not sensitive to viscous effects, for 
example the flow past a square body at Reynolds numbers between 10 4 and 
10 7 . KPD 1 has been successfully usti for the "Vortex Flowmeter” study 
with Dr. Couet [40] . This configuration involves several interacting bluff 
bodies. 

KPD2 is directly based on KPDl; it treats the boundary layer, from the 
attachment point to the first separation point, with an integral method and 
treats the rest of the domain with the Vortex Method. This removes the 
problems with premature separation experienced with KPD 1. The integral 
method is designed for boundary layers imbedded in an irrotational flow; 
moreover, it exhibits a singularity at the separation point. This is why 
it cannot be applied beyond the separation point. The boundary layer is 
also considered as quasi-steady. KPD2 is suited to problems with a single 
streamlined body, and can capture the major viscous and turbulent effects to 
which K PD 1 is insensitive. It has been used mostly for airfoil flows, including 
dynamic stall [39] and the tilt-rotor configuration (work to be published in 
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cooperation with W. J. McCroskey). 

KPD3 is quite different from KPDl or KPD2) it is the latest version of 
the program and possesses the most potential. It treats the viscous region 
along the wall with a truly unsteady implicit finite difference boundary layer 
method, in a manner that is valid even inside a vortical outer flow (like the 
wake of the body itself or the wake of another body). The boundary layer 
solver is not classical. It allows for separation and reattachment of vorticity; 
intuitive arguments are invoked to couple it as strongly as possible with the 
outer solution. K PD3 has been validated on the circular cylinder at moderate 
and high Reynolds numbers, using the Baldwin-Lomax algebraic turbulence 
model. It can treat several bodies without special precautions, and in general 
is more accurate and provides more information than KPDl or KPD2. It is 
not quite as robust in its present version; in particular it can have difficulties 
near sharp edges. 
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H ANALYTICAL CONSIDERATIONS. 


1) Vorticity formulation. 

The incompressible Navier-Stokes equations are formulated in terms of the 
vorticity, and the resulting system of equations will be solved numerically. It 
will be shu*n that the initial-value problem used for the vorticity is mathe- 
matically equivalent to i he i< ual initial- value problem used for (u, v, p). Since 
the "(u, v,f," system is well posed, the "w” system will also be well posed. 
The vorticity formulation is considered to be more efficient numerically. 

The first subsection will introduce the necessary definitions and present the 
formal proof of equivalence of the two systems. The procedure follows closely 
the work of J. Wu, described in [14]. The second subsection will contain 
some comments about the aspects of the procedure that do not follow the 
traditional train of thought and sometimes become misunderstood. 

a) Definitions and proof of equivalence. 

The domain is the ( 2 , y) plane. It contains M solid regions called S m ; each 
Srn is an open and bounded domain with a boundary dS m ■ Let 3 be the 
union of the S m 's and F be the fluid domain. Thus the plane is partitioned 
into the two open domains S and F and the boundary dS. In general, the 
solids move and therefore S and F depend on the time, t. 

The momentum equation, Eq. 2, contains a parabolic diffusion term, i/AU, 
and therefore the function U is expected to be smooth. U is considered as 
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being at least of class C3 in F, that is three times continuously differentiable, 
at all times t > 0. The pressure is at least of class C\. 

The velocity of a point (x, y) belonging to the solid S m is given by the 
function U m : 

U m (i,y) = U m „ + n m (“ ! ') (8) 

where U^o and tt m are known functions of time. In this study the motion 
of the solid bodies will always be prescribed, but the theory would not be 
different if it were known from the solution of a dynamic equation (for 
example a solid with elastic restraint). 

The complete system of equations governing u, v and p is the following: 


(incompressibility, Eq. 1) V.U = 0 in F (9.1) 

(momentum, Eq. 2) = —Vp + 1 / AU in F (9.2) 

(at wall) U(2, y) = V m (x, y) on dS m (9.3) 

(far field) U f Joo for ||r|J 00 (9.4) 


(initially irrotational) At f = 0 V x U = 0 in F (9.5) 

(no initial circulation) At t — 0 f V.ds = 0 (9.6) 

JCm 

where C m is a contour that encloses S m and dS m . A more accurate definition 
of the contour within F is not necessary, because the velocity field is irrota- 
tional in F at time zero (Eq. 9.5), so that the line integral does not depend 
on the contour. 

Let us turn our attention to the vorticity formulation. The vorticity u is 
defined by Eq. 3. Since U L considered as being of class C3, u is considered 
to be of class C 2 in F for t ^ 0. In the exact solution the vorticity is known 
to decay exponentially at large distances from the body, provided that it 
did initially, at time zero [14]. All the flows considered here will be started 
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from rest, with zero vorticity in F; therefore exponential decay of u can be 
assumed. As a consequence, all generalized integrals involving the vorticity 
over the infinite region are absolutely convergent and have the same behavior 
as if the vorticity had a finite support. 

For the vorticity formulation it is convenient to formally extend the velocity 
field to cover the whole plane; inside F it is the fluid velocity and inside S m 
it is the velocity U m of the solid material. Naturally, the same dynamical 
equations do not apply in F and in S (in particular the pressure will not 
be extended to 5), but this does not affect the kinematics of the extended 
velocity field. The reason for extending the various fields into S is that 
this will allow the use of Green’s functions without any boundary terms or 
"images” for the solution of the Cauchy-Riemann equations. 

Similarly, an extended vorticity field is defined by applying the definition, 
Eq. 3, both in F and in S. Inside S the velocity (given by Eq. 8) and the 
vorticity are both of Class Cqq. 

We can *;ow introduce the system of equations that will govern the vorticity: 


(vorticity conservation law, Eq. 4) 

Du 

~ = i/A u in F 

Dt 

(10.1) 

(vorticity inside solid) 

u = 2 O m in S m 

(10.2) 


(Biot-Savart) 



-f-oo -hoc 

— f f ~ 

2ir J J \x — x‘){x— x'f + {y — y') 2 

— 00 — OO 


(10.3) 


(at wall) U.n = U m .n on oS m (n : normal vector to d S m ) (10.4) 

(additional condition) f — —2 Am——- ( A m : area of S m ) 

J on at 

9S m 

(10.5) 

(initially irrotational) at t = 0 u = 0 in F (10.6) 
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(no initial circulation) at t = 0 


// 


ui(x,y)dxdy = 0 


(10.7) 


D m 

where the domain D m encloses S m and dS m like the contour C m in Eq. (9.6). 
The main result of this chapter is that the systems 9 and 10 are equivalent. 


First let’s prove that the system 9 implies the system 10: 

• Equation 10.1 has already been derived as the curl of Eq. 9.2. 

• Equation 10.2 is obtained by taking the curl of Eq. 8. 

• Equation 10.3 is the Green’s function integral giving the solution U of 
the Cauchy Riemann system formed by Eq. 9.1 and Eq. 3, subject to the 
boundary conditions, Eq. 9.4. 

• Equation 10.4 is a consequence of Eq. 9.3. 

• To derive Eq. 10.5 we use Eq. 9.2 and the identity: 


AU = V(V.U) — V x w 


in addition V.U = 0 from Eq. 9.1. We then write Eq. 9.2 on dS m and take 
its dot product with the tangent unit vector -. 


DU , dp . dcu 

— — .t v— 

Dt ds dn 


(12) 


Since the particles, locally, adhere to the wall their acceleration is the same 
as the acceleration of the wall: DV/Dt — DV m /Dt on dS m . Thus we have: 


DUm , _ _dp . du 
Dt ' ds V dn 


(13) 


We then integrate Eq. 13 along the closed contour dS m . The acceleration 
derived from Eq. 8 is integrated analytically, and the pressure term cancels. 
The final result is Eq. 10.5. 


• Equation 10.6 is a consequence of Eq. 3 and Eq. 9.5. 
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• Finally, Eq. 10.7 is a consequence of Eq. 3, applied in F and S m , and 
Eq. 9.6. 


Let us now prove that the system 10, in return, implies the system 9: 

• Equation 9.1 is automatically satisfied when the velocity field U is 
generated by the Biot-Savart law, Eq. 10.3. 

• Equation 9.4 is also automatically satisfied due to Eq. 10.3 and the fact 
that cl/ decays exponentially away from the origin. 

• To prove Eq. 9.3 it is convenient to introduce a stream function. The 
velocity U given by Eq. 10.3 is divergence-free and a stream function ip can 
be associated with it and given by: 


dtp 

dx 


= v 



(14) 


The solid body velocity field given by Eq. 8 is also divergence-free; a stream 
function tp m can be associated with it and defined over S m by the same 
formula as Eq. 14. In both cases the stream function is defined except for an 
arbitrary additive constant. Then Eq. 10.4 can be rewritten: 


dip _ dip m 
ds ds 


along 


dS m 


(15) 


Thus ip — tp m is constant along dS m . In addition, as a consequence of Eq. 3 
and Eq. 14 the following Poisson’s equation applies: 


A ip = ui 


(16) 


Now the (scalar) curl of U m is 2Q m , and u is also equal to 2fi m in S m , from 
Eq. 10.2. Therefore ip and tp m satisfy the same Poisson’s equation, and: 

A(V> — tp m ) — 0 in S m (17) 
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The function ip — ip m satisfies Laplace’s equation in S m , which is bounded. 
It is also constant on dS m (Eq. 15), which represents a Dirichlet boundary 
condition. This is a well posed problem and the unique solution for ip— ip m is 
a constant over S m . Therefore its derivatives are zero, which can be written 
as: 

U(z, y) = U m (x, y) in S m (IS) 


The velocity field is equal to the solid body velocity inside the solid. 
Furthermore since u> is considered to be of class C 2 it is bounded (for t > 0) 
and the velocity field U generated by Eq. 10.3 is continuous; so if it is equal 
to U m inside S m it is also equal to U m on dS m , (the solid body is assumed 
to have a finite thickness) and Eq. 9.3 follows. 

• To prove Eq. 9.2 we have to produce a pressure field. Let us consider 
the quantity: 

52-i/AU (19) 

Dt v 

with U given by Eq. 10.3. If we take the curl of Eq. 19 we get (s'uce U is 

divergence-free): 

< 20 > 

which is zero from Eq. 10.1. Therefore the quantity defined by Eq. 19 is the 
gradient of a function p: 


i/AU = -Vp 

Dt y 


(21) 


We now write Eq. 21 on dS m , rewrite the viscous term as in Eq. 13, use the 
same argument for DU/Dt and take the dot product with the tangent unit 
vector to obtain: 


/ = + * / Tn dS (22) 
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The right-hand side is zero from Eq. 10.5. This means that p is single-valued 
around each solid S m . Therefore p is the pressure (defined except for an 
additive constant) and Eq. 9.2 is satisfied. 

• Equation 9.5 is a consequence of Eq. 10.6 and Eq. 3. 

• Finally, Eq. 9.6 is a consequence of Eq. 10.7 and Eq. 3, app’ied in F 
and S m . 

This completes the proof. 

b) Comments. 

The first comment qualifies out assertion that "the system is well-posed". 
The far fieldcondition, Eq. 9.4, is imprecise in the sense that it does not 
specify how fast the difference (U — U,*,) tends to zero as ||r|| tends to oo. 
How strong this decay should be to provide a well-posed system with the 
Navier-Stokes equations has not been rigorously established. The common 
practice one follows, when confronted with this question, is to as ume a 
behavior in the far field that is as regular as possible. If we assume that 
the velocity can be expanded in negative powers of ||r|| and that the tlow 
is effectively irrotaticnal in the far field (the vorticity decays expouentially) 
then the terms of order ||r||~ 1 are a source term and a vortex term. The 
source term must be zero for mass to be conserved. The vortex term gives 
the circulation around a large contour. This circulation must be a constant, 
from Kelvin’s theorem (the viscous term t'AU has been written — rV x u ; and 
therefore decays exponentially, if u does). If a steady lifting flow is sought 
the circulation will not be zero. In our case the value of the circulation 
does not matter much since the flow is viscous and unsteady, and thus will 
wash away any excess circulation. We shall assume zero circulation; therefore 
the velocity disturbances decay like ||r|| — 2 . With this decay specified, the 
Cauchy-Riemann system has a unique solution, given by Eq. 10.3. Rigorously 
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this should not be considered as a pure boundary condition: it would be over- 
specified. It contains the boundary condition and an assumption about the 
far-fleld behavior of the solution. 

The initial conditions were derived in the same spirit. They simulate an 
impulsive start of the flow. After such a start the vorticity will be confined 
to dS and be infinite in magnitude, with a finite jump of velocity across <95. 
This is simply a potential flow problem and it is well known that in such a 
case the circulation around each solid is arbitrary. It was set to zero. 

The second comment concerns the boundary condition at the wall. If 
we examine the system 10 and especially Eq. 10.4, it seems that only the 
continuity condition (zero velocity normal to dS m ) is applied and that the 
no-slip condition (zero velocity parallel to <95 m ) has been lost. However it was 
shown that Eq. 9.3, which includes no-slip, was satisfied. This paradox is 
clarified by noting, first, that the velocity fields produced by the formula 10.3 
are not arbitrary (they are divergence-free and have the required vorticity 
2n m in 5 m ), and second, that what we have shown is that the global normal 
velocity condition (Eq. 10.4 applied all along <95 m ) implies the global no-slip 
condition. Naturally, the local normal velocity condition does not imply the 
no-slip condition. 

Except for Chorin’s first paper [29] and the recent paper by Lewis [34], all 
papers employing vortices imposed the boundary condition, Eq. 9.3, in two 
stages. First, they included image vortices in the Biot-Savart law, Eq. 10.3, 
to secure the normal velocity condition, and then, they introduced vortices 
to secure the no-slip condition. Here the complete boundary condition is 
obtained in one step by introducing vortices to satisfy Eq. 10.4, which 
becomes an integral equation for co if U is replaced by Eq. 10.3. This 
is much more efficient since Eq. 10.4 and 10.3 can be written directly in 
the physical plane, whereas the image vortices could only be used after a 
conformal transformation had converted the body into a circle. This made 
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the Vortex Method awkward and limited the simulations to the few cases for 
which the mapping and its derivatives are . nown: ellipses, Joukovsky airfoils, 
etc. 

The third comment is about the conservation of circulation. We mentioned 
previously Kelvin’s theorem concerning the circulation around a large con- 
tour, T. Equation 10.5 has an interesting consequence which extends Kelvin's 
theorem. This result is due to Wu [14]. T is equal to 

+oo -f-oo 

’-II u i{x,y)dxdy (23) 

— 00 — oo 

(the integral includes the vorticity that is ins id* S). To evaluate dT/dt it is 
convenient to use Lagrangian coordinates, because in Lagrangian coordinates 
F and 5 do not depend on t, so that points do not switch from F to S as the 
solids move. Let (c, b) be Lagrangian coordinates which coincide with ( x,y ) 
at the time considered. The Jacobian of t i > mapping (o,6) -*• (X,y) is 1 at 
any time since the flow is incompressible. Therefore T i s also equal to: 

-foo -foo 

'-/ / u(ji,b)dadb (24) 

— 00 — OO 


we can integrate in either set of coordinates. Then dT/dt is: 

ft = / / ^ a ^ dadb < 25) 

— OO —00 

This is Reynolds’ transport theorem. 

We can now revert to the (x,y) coordinates to evaluate Eq. 25. In S m 
Du/Dt is 2dCl m /dt and integrates to 2A m dTt m ldt. In F, Duj/Dt is v and 
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is easily integrated by parts, to yield: 



There is no contribution from the far-fleld, because of the exponential decay of 
u>. In addition, for each m the expression between brackets is zero, according 
to Eq. 10.5. The integral along dS m that appears in Eqs. 10.5 and 26 can 
be interpreted as the total production of fluid vorticity along dS m . Clearly 
each solid, while changing its internal circulation at the rate 2A rn dll rn /dt . , 
releases an opposite amount of vorticity into F. If there is only one solid this 
is equivalent to Kelvin’s theorem, which states that dT/dt is zero. If there 
are several solid bodies Eq. 26 is a stronger result, since each solid separately 
contributes zero to the circulation. 

Another point of interest is the way the pressure is computed. Computing 
the pressure is not necessary in order to solve the vorticity equation, but it 
is an excellent way to monitor the simulation. The common wj v to interpret 
boundary layer behavior is in terms of the pressure gradient along the wall. 
Formally, the pressure is given by Eq. 21: however this equation would be 
hard to use numerically with the Vortex Method. On the other hand Eq. 13 
gives the wall pressure gradient as a function of vdu/dn, and we have seen 
that vdu/dn is the rate of creation of vorticity at the wall; this quantity 
is well defined in the Vortex Method and will allow the wall pressure to be 
computed accurately, even inside the wake. Using Bernoulli’s equation in the 
wake would be incorrect since tne flow there is vortical. 

A detail remains: Eq. 13 only yields the pressure gradient; thus the pressure 
is known except for an additive constant. If one wishes to determine this 
constant and the body is in contact with the irrotational region, it is possible 
to apply Bernoulli’s theorem from infinity upstream to a point on the attached 
part of the boundary layer. In practice it is convenient to use the front 
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stagnation point. 

Another useful boundary layer monitor is the wall shear stiess. It is equal 
to: 

t — vu on dS m (27) 

2) Approximations for the outer and inner regions. 

Matching procedure. 


a) Motivation. 

This section describes the approximations that are made aDd how certain 
considerations allow us to simplify the equations, by omitting terms that are 
known to be small or information that is not important. 

Moot of this section applies to all three versions of the program; when they 
differ, the description will apply only to KPDS. The theory implemented : n 
KPDl and KPD2 and their numerical aspects will be described in Appendix 
A. 

The most important, and the most delicate, approximation is naturally the 
neglect of viscosity. The coefficient of viscosity is small, but it multiplies the 
highest derivative, and the perturbation problem is said to be singular [42] . 
The inviscid problem and the viscous problem have very different characters; 
in particular they do not require the same number of boundary conditions. 
Regions exist in the flow where the velocity gradients are so large that the 
viscous term is as large as the inertia term. This vi r oous term can change the 
locaf \alue of the vorticity by an amount of order 1, meaning that it does not 
tend to zero while the coefficient of viscosity docs. Therefore the flow with 
small viscosity cannot be treated as slightly different from the inviscid flow 
in the usual sense, and a straightforward attempt to expand the solution as 
a power series in v would fail. 
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The justification for omitting the viscosity is the following. The effect of 
viscosity will be to diffuse the vorticity over very short distances, without 
creating or destroying any. (The viscous term in Eq. 4 is the divergence 
of i/Vuj and uVu is interpreted as a flux of vorticity. It is not a source 
term.) Let L be the length scale associated with the body, Uoo the freestreani 
velocity and 1 / the kinematic viscosity. The non-dimensional number LVoofu 
is the Reynolds number, and is large in all cases under consideration (over 
10 4 ). The length scale associated with the viscous diffusion is \ful where t 
is the ” age” of the vorticity. Let us consider some vorticity which is ” born” 
at the solid boundary and in a time L/ Uoo is transported into the wake, to a 
distance L from the solid. The viscous scale becomes \Z r L'L/ij 00 and the ratio 
of this scale to L is y/vpJJ^, or Re~ 1 / 2 , and thus is small. Integrals like 
the one in equation 10.3 and in general the flow close to the solid boundary 
will not be sensitive to the displacement of the vorticity over such a small 
distance. Since predicting the stresses on the solid is the ultimate objective 
of the study, omitting detailed information about the vorticity diffusion in 
the wake is minor as long as the transport is correct. 

However the vorticity is "produced” at the solid boundary [43] and its 
subsequent transport is very sensitive to its initial life, near the wall, during 
which the scales are small and the viscous term important. It is the convection 
with the fluid that carries the vorticity into the large structures of the wake, 
but the velocity is zero at the wall and only the viscosity can make the 
vorticity penetrate into the stream at all. Therefore the "justification” we 
just reviewed breaks down in the wall region. 

This motivates the procedure, illustrated in Fig. 5, of coupling an inviscid 
outer flow and a viscous boundary ! \yer flow. This procedure is common when 
the outer flow is not only treated as inviscid but also as irrotational. Here, 
the outer flow will be vortical. The effort will be worthwhile if an efficient 
solver is available for the simplified equations in each region. 
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The Vortex Method is efficient in the outer flow; it treats the transport 
term accurately and provides the necessary resolution in the wake. It does not 
cause any problem at large distances from the body. Its weakness in treating 
the detailed viscous features will not disturb the large scale structures which 
dominate in the wake. The implicit finite difference method is very good 
at treating thin viscous flows. For such a small and logically rectangular 
domain it is also very fast. Both methods are available and well tested. The 
new element that is needed is a procedure that makes the two regions interact 
through the boundary conditions at the interface. 

In a previous investigation Shestakov also coupled the Vortex Method to 
an Eulerian method [44] ; however he used the Vortex Method in the wall 
region, and the Eulerian mu, hod away from it! Even though the conditions 
were slightly different (he treated an internal flow) our interpretation and 
Shestakov’s appear to be totally opposite. His results appear reasonable, but 
it is not clear how much his flow depended on the wall region, or how much 
benefit he derived from using the Vo iex Method near the wall. 

The two approximate systems of equations will be described separately, 
followed by a discussion of the conditions at the interface. 


b) Outer flow. 

In the outer region, the viscous term in Eq. i is dropped, only the transport 
term is retained. The approximate equation is: 


Du 

~Dt 


= 0 


(28) 


The material derivative of the vorticity, or equivalently its time derivative 
in Lagrangian coordinates, is zero; this is what makes a Lagrangian method 
attractive. 

The vorticity is zero at large distances (the system is always started with 
the fluid at rest) and no boundary condition is needed in the far field for Eq. 
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28. The far field behavior of the velocity is essentially the same as in the 
exact formulation, in which the vorticity decayed exponentially. 

The proper boundary condition for the hyperbolic Eq. 28 at the interface 
with another domain depends on the sense of the velocity; u i itself is the 
characteristic variable and the velocity U gives the characteristic direction. 
If this velocity is into the other domain (outflow), no condition should be 
applied; if it L into this domain, the value of the vorticity, or equivalently 
the flux of vorticity, should be prescribed. Information travels in the same 
sense as the particles, and this is realized very simply with a Lagrangian 
method: an outflow boundary absorbs particles and information, an inflow 
boundary generates new particles which carry information. 

c) Inner flow. 

In the boundary layer the viscous terms are retained, but the thinness of the 
layer renders some terms negligible. Curvature effects will not be included. 
This is legitimate for shapes like a circular cylinder; for airfoils, it might be 
necessary to account for curvature near the trailing edg , or to round it off 
so as to increase its radius of curvature. 

Let s and n be the coordinates along the wail and normal to it respectively, 
and u and v be the velocity components in the s and n directions respectively. 
The scale in the n direction being much smaller than that in the s direction 
allows Eqs. 1,3 and 4 to be approximated by: 



du , dt) 

(29) 


=■ 0 

ds ^ dn 


du 

U dn 

(30) 

du 

d 2 ui 

(31) 

dt 

+ U.Vuj = r— 


Equation 1 has simply been reformulated in terms of (s,n,u,v), without 
approximation, to yield Eq. 29. The definition of the vorticity, Eq. 3, has 
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been simplified by dropping the dvjds term, giving Eq. 30, and the viscous 
term in the s direction has been dropped from Eq. 4 to give Eq. 31. 

These differential equations are the same as the conventional time- 
dependent boundary layer equations, but the boundary conditions employed 
■will be different. The inner region extends around the whole body, and the 
boundary conditions in the s direction are periodic. The equation is advanced 
in time, not by ” marching” along the boundary layer in the direction of the 
local velocity. Thus the type and stability of the equation are not affected 
when this velocity changes sign, for instance at separation. The other major 
difference is that, whereas conventional boundary layers are imbedded in an 
irrotational outer flow, this one is not. In particular, there is no Bernoulli 
relation linking the outer velocity to the pressure gradient. Also, the vorticity 
does not necessarily tend to zero at the outer edge of the inner region, and 
the boundary condition at this edge must allow a transfer of vorticity to or 
from the outer flow. 

The edge of the inner region is at n — 6, where <5 is a parameter. 6 
should be small enough for the boundary layer assumptions to be valid; on 
the other hand 5 should be large enough for the physical viscous region to 
be contained in the computational region. Naturally, the "viscous region’' 
cannot be precisely defined; however, if the inner solution reveals strong 
gradients confined to the vicinity of the wall and a quieter region elsewhere, 
6 is probably large enough. Another test is tc compute the various physical 
thicknesses of the boundary layer (displacement, momentum, etc ) and to 
compare them to 6. Along the attached region, the boundary layer is well 
within 6; after separation almost all the vorticity is in the outer region and 
thcie is no boundary layer in the usual sense. Examples will be given to 
illustrate how 6 is chosen. 

In the boundary layer the velocity is obtained by integrating Eqs. 29 and 
30 in the n direction. Both components of the velocity are zero at the wall. 


37 



This provides a well posed system with Eqs. 29 and 30 since these are first 
order. 

Equation 31 is of second order in the n direction and thus requires two 
boundary conditions. 

Instead of a boundary condition at the wall, the integral equations, Eqs. 
10.3, 10.4 and Eq. 10.5, are used. This is natural since we have shown that in 
the exact formulation Eqs. 10.3, 10.4 and 10.5 regulate the flux of vorticity 
from the wall. 

The other condition regulates the flux of vorticity through the interface. 

d) Interface conditions. 

Both the velocity and the vorticity field should be matched at the interface. 
The matching of the velocities is done in the same spirit as in classical 
boundary layer theory. The component of velocity parallel to the wall, u, 
will always be matched since it is of order 1. The normal component v is 
small, of order S, at n = S and at the lowest level of approximation it is 
neglected; we shall adopt the next level of approximation and match the v 
components as well(still assuming that the inner region is thin). 

As for the vorticity, since Eq. 28 is first order and Eq. 31 is second order, 
they cannot be matched without making an additional approximation. 

The two domains exchange vorticity through the interface. Since Eq. 28 
is applied down to the interface, it is consistent to derive the approximate 
interface condition in the same spirit. Thus it is assumed that the interface 
is far enough from the wall for the viscous term to be dominated by the 
convection term, and the transfer of vorticity is taken as vu and imposed by 
the upstream region. The boundary condition has thus dropped to ' he level 
of the inviscid approximation. For Eq. 31, this means that the viscous term 
neglected at the outer edge of the boundary layer. 
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HI) NUMERICAL IMPLEMENTATION. 


Discrete approximations to the continuous equations are derived as a basis 
for the following discussion of the numerical method. The algorithm used 
for each region are described first, and then the coupling procedure is intro- 
duced. Each method converges in its domain as the scale of the discretization 
is reduced in space and time. However, the complete algorithm should not 
be expected to converge to the Navier-Stokes solution since the errors intro- 
duced by the inviscid and the boundary layer approximations remain finite. 
Convergence at a given Reynolds number could only occur if the order of the 
boundary layer approximation (among other things) was increased in parallel 
with the numerical refinement. 

This chapter applies to the KPD3 program; KPD1 and especially KPD2 
use a different logic which will be described in Appendix A. 

1 ) Outer flow. 

a) Extent of the outer region and discretization in space. 

The outer region covers the whole ( 2 , y) plane except the solid and a narrow 
band of thickness 6 around it. It extends to infinity and no grid is involved. 

The vorticity field is described as the sum of a large number N v of mobile 
functions of small support, referred to as "vortices” . Each vortex is defined by 
the position r* = (a:,, y») of its center, its circulation I\ which is the integral 
of the vorticity it carries, and the shape 7 of the distribution of the vorticity 
around the center (see Eq. 5). This distribution is in general bell-shaped; this 
is the "vortex blob” method. 
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Ail individual vortex does not live for the full duration of the computation 
New vortices are created at the interface, where they enter the outer flow. 
Vortices can also be absorbed by the wall region and thus removed from the 
computation. Finally, vortices are allowed to merge when certain conditions 
are fulfilled. The vortices are independent entities; they form a "jungle” 
except maybe just after separation, where the free shear layer has not yet 
undergone instability and broken down into circular eddies. 

In this study, the shape of the blob, 7 is taken to be axisymmetric, and is 
the same for all vortices at all times. The whole blob moves at the velocity 
of its center. Clearly, no diffusion of the blob takes place, and the straining 
of the blob by the velocity gradients is also neglected. This straining is the 
source of the spatial error, as analyzed in [9] and [19]. 

Simple cores, defined by algebraic functions, were used. Being everywhere 
positive, these cores are expected to yield second order convergence in an in- 
viscid problem [19]. The superiority of the more elaborate cores (which should 
yield higher order convergence) has not been clearly demonstrated [45] ; there- 
fore the simplest possible approach was chosen. The cores chosen also require 
less computing time. The computation of the interactions is the most time- 
consuming part of the program and it might be advantageous to have many 
" inexpensive” vortices rather than a smaller number ~>f ” sophisticated” ones. 
Finally, it is very likely that the main source of error is not in the treatment 
of the inviscid transport of vorticity, but in the neglect of the small scale 
turbulence, and even more in the interaction with the walls. In the wall 
region the solution is not smooth at the scale of the vortex core radius, and 
the rate of convergence of the method becomes less relevant. 

Two cores were used and are defined by 


Core 1 
Core 2 


7 (r) 


-( 2 # 

lo, 


if r < cr; 
if > cr. 


7(f) 2rr(r 2 -f j 2 ) 2 


(32.1) 

(32.2) 
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a is the "core radius". These functions are plotted in Fig. 6, as well as 
the corresponding velocity and stream function distributions, with the point 
vortex as a reference. In Core 1 the vorticity is continuous and the velocity 
continuously differentiable, while Core 2 is infinitely differentiable. Core 1 has 
a finite support and a vortex is not allowed within a distance a of the wall, 
so that the vorticity is entirely outside the solid. Core 2 does not have this 
property: its support is i nfini te (the vorticity decays like r ~ 4 ) and penetrates 
the solid. Although this does not appear to be very natural, the differences 
in the results were negligible. 

Core 1 was used in some versions of KPDl and KPD2, in particular for 
the case of dynamic stall (on a CDC 7600). It was then decided to switch 
to Core 2 for the CRAY version of the program, because Core 1 involves an 
"IF” test which inhibits the vectorization of the loop. 

b) Computation of the velocity. 

The velocity field must be computed in order to solve Eq. 28 for w. If we 
introduce the value of u> from Eq. K into the Biot-Savart law, Eq. 10.3, the 
velocity induced at a point r by the vortices is given by: 



with 7) defined by: 

^JLUl — r7(r)and r\ sa r“ 2 for rlarge. (34) 

or 

The formula is greatly simplified by the fact that the blobs are axisym- 
metric. If point vortices were used r] would be equal to r“ 2 . With vortex 
blobs rj is regular near zero, and the velocity field is smooth. 

The uniform freestream velocity also contributes to the velocity field, as 
well as the velocity induced by the inner flow vorticity, which is not included 
in the blobs. Since the inner region is a thin shear layer, even compared to the 
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core diameters of the vortices, it can be represented to a good approximation 
by a vortex sheet of zero thickness. The strength of the vortex sheet will 
be designated by U e ; it is the circulation per unit length of the sheet, and 
also the jump of tangential velocity across the sheet. The sheet is made up 
of segments, each segment covering the interval between two wall points (see 
Fig. 7). The strength of the vortex sheet is assumed piecewise linear. The 
velocity field of such a segment of vortex sheet can be expressed analytically. 
In complex variables it is 


U(*)« 


i ( z 2 - z\) 


2 ir \z 2 — z i| 


V* ~ U e2 + 


Ueljz — *2)— UMZ — fQ 

(z 2 — Z\ ) 


log 


( - — * 2 Y 

\z - Zl). 


(35) 

where z\ and z 2 are the two ends and t7 el and U e2 are the strength at 
each end and the overbar denotes the complex conjugate. This field jumps 
across the sheet but is smooth on each side: this is why segments are used 
instead of circular vortices to represent the inner flow vorticity. The velocity 
field however has a logarithmic singularity at the junction of two adjoining 
segments unless they have the same slope. Therefore it is desirable to keep 
this slope as smooth as possible. 

Thus the velocity of each vortex is the sum of the freestream velocity Uqo, 
N w terms of the type given by Eq. 35 for the N w wall intervals, and N v 
terms of the type given by Eq. 33 for the JV„ vortex blobs. This is the 
discrete analog of Eq. 10.3. 

The computation of the interactions has to be performed at each time step 
and this is the most time-consuming part of the program. N V {N V + N w ) 
interactions have to be computed, and each of the N V N W interactions with 
the wall segments involves a complex logarithm. Fortunately, the simplicity 
of the data base makes vectorization easy, provided that the function rj does 
not involve "IF" tests or any non-simple function. Even then, it is worthwhile 
to apply analytical tools to reduce this cost. 

The high cost of implementing the Biot-Savart law, Eq. 33, comes from 
the fact that each vortex interacts with vortices in the whole domain, with 
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distant vortices as well as with its neighbours. On the other hand, the velocity 
field induced by a distant cluster of vortices does not depend much on its 
detailed shape, and this should be taken into account by the program. The 
velocity induced at a large distance R by a cluster of diameter / has a Taylor 
expansion in terms of l/R. The first terms of this expansion are a vortex 
term, a dipole term, and so on. 

To implement this in a controlled way, with a known and bounded error, 
the clusters are first given a precise definition. The plane is divided into a 
number n of identical square cells of side /, surrounded by an external cell 
which is treated separately (Fig. 8). For the scheme to achieve its purpose, 
each cell should contain more than a few vortices; so n should be much smaller 
than N v . Each time the interactions are to be computed, the vortices that 
are in the same cell are linked, logically. Their distance to the center of the 
cell is smaller than l/>/2. 

It is convenient to use complex notation here. The function that is ex- 
panded is (z, — Zj)~ 1 , where z, is complex for (x, ,y % ). Let z, be in the K th 
cell, with center Zk, and z ; in the L th cell (see Fig. 9). Thus (z< — z ; ) — 1 is 
expanded in the vicinity of (Zk — Zl)~ l - The function z~ l has a rapidly 
converging Taylor expansion, and the error can be bounded as a function of 
1/\Zk — Zi\ and of the number of terms that are retained. This number of 
terms is chosen to make the error as uniform as possible. If the two cells 
are far from each other, compared to /, only the first term of the expansion 
will be kept. If they are not very far, up to four terms will be included. If 
they are neighbours, the Taylor expansion does not apply; in that case, the 
Interactions are computed vortex by vortex. If either vortex is in the external 
cell, the Taylor expansion is not used either, since the external cell is infinite 
in extent and has no "center". 

In the final version of the program, enough terms were taken to ensure 
a maximum relative error of l?o in each interaction. The actual error was 
computed in a test case by also computing the velocities without using Taylor 
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expansions. The maximum difference that was observed was about 0.59c of 
Uoo, and tu? " Lz average was less than 0.1%. This level of error is very 
moderate, >n comparison to the other possible sources of error. 

With a proper choice of / and n, computing the interactions this way 
instead of applying Eqs. 33 and 35 directly can save 65% of the time on 
a serial computer like the CDC 7600. Typical values are N v = 1000, n = 
100, / = .5 with a oody of size 2. Vectorizing the Taylor expansion for the 
CRAY was possible, but resulted in a program that was more complex and 
less "smooth" logically, with shorter loops. As a result, KPD1 and KPD2 
run faster without the Taylor expansions; KPDZ still runs faster with them, 
because using them saves the time of evaluation of most of the complex 
logarithms in Eq. 35. 

c) Time integration. 

The system of ordinary differential equations, Eqs. 7 and 24, is integrated 
by the Adams- Bashforth second order method. The velocity of each vortex 
is computed at uniform time intervals and the positions updated according 
to the formula: 

r.(t + AI) = r.(f) + At(jU.(f)- Iu,(i _ At)) (36) 

At is the time step and « .e accuracy in terms of At is of second order [46] . 

As with any multistep method, the first step of integration mu it be treated 
differently because U*(t — At) is not available. Thus the first step in the life 
of each vortex is handled by the explicit Euler scheme: 

r,<t + At) = r,(t)-f AtU,(t) (37) 

The Adams-Bashforth scheme was chosen because it is second order accurate, 
while requiring only one evaluation of the derivative per step. It is weakly 
unstable when applied to linear equations, but the non-linearity of equations 
7 and 24 actually stabilizes the integration and no stability problem has been 
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encountered (see subsection e)). The need to store two levels of the velocity 
values is not a problem since the Vortex Method involves only a very moderate 
number of variables. 

It should be noted, however, that since the inner region solution is only 
first order accurate in time, the overall accuracy is of first order at best. 
The Adams-Bashforth schema is used mainly to gain quantitative accuracy 
(over the first order Euler scheme) and especially reduce the scattering of the 
vortices (see subsection e) ). 

d) Vortex merging device. 

The boundary layer releases a significant number of vortices near the wall 
at each time step: typically 100 new vortices, compared to a total number of 
1000. (Naturally, these luO vortices do not carry 10% of the vorticity; they 
are numerous but weak. Typically, in one time step the new vortices of one 
sign might add up to a circulation of 0.03, while each one of main "Kerman 
Street" eddies carries a circulation of the order of 10). This continuous 
addition of new variables should be balanced by the suppression of some of the 
old variables at approximately the same rate; this is done by merging pairs 
of vortices into one when appropriate conditions are fulfilled. As a result, 
the vortices are dense near the body, where fine resolution is desirable, and 
become progressively sparser away from it. 

Deffenbaugh and Marshall introduced a merging method l it did not make 
all the details available [28]. R. Rogallo (personal communication) used a 
device which was very similar in spirit to the one used here; however the 
error estimate was different. 

If nothing was done to keep the number of vortices under control the 
program would only be able to compute flows of relatively short duration 
before the number of voriices and the associated computing cost would be- 
come excessive. This would be acceptable for some applications (slender body 
"2D-3D" analogy, for instance) but not for the ones considered in this study. 
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Furthermore, if all the vortices are kept, there is a strong incentive to 
create fewer at each step; a* a result, the wall boundary condition is more 
loosely satisfied: typically only 20 discrete equations are retained [29]. It 
is preferable to have a much greater resolution near the body, typically 200 
discrete equations, and then let these many small vortices progressively merge 
into larger ones. Furthermore the description of the flow is more homogeneous 
in time and .’an actually reach an asymptotic state. 

The procedure is the following. Only the merging of two vortices into one 
is attempted, at each time step. If we consider two vortices of circulation r t 
and r 2 and position z\ and z 2 , the velocity field they create before merging 
is, in complex notation: 


U(*)« 



£i 

(z - z x ) 


+ 



{z — z 2 ) 


The field they create after merging is: 


uV) = 


■ r 

2 *(*- Z) 


(38) 


(39) 


where T and Z are the circulation and position, recpectively, of the new vortex. 
The first few terms of the expansion of the difference U(£)— U'( 2 ) at large 
distance z are (the complex conjugate of): 


• g r - r, - r a ) 

2 


0V, + r 2 j 2 -rz) (rz a -r,^-r 2 4) 

..O I Q 


) 


+o(kr*) 


The two leading terms can be removed by taking: 


(40) 


r = r, + r 2 

r*i ~ i -f 


(41.1) 

(41.2) 
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This means that the new vortex is given the sum of the circulations of the 
old ones and placed at their centroid (Fig. 10). It is worth noting that 
two vortices, if they are isolated, orbit precisely around this centroid, which 
is itself stationary. Thus by merging we replace two vortices, which would 
move on two concentric circles, by one stationary vortex at the center of 
these circles. This way of merging also preserves the total circulation, and 
the first moment of vorticity, which is equal to the impulse of the flow (this 
results from an integration by parts [14]). It should also be noted that if the 
two vortices have opposite signs the centroid is aligned with them but not 
between them (Fig. lO.b). 

The third term cannot be removed within this framework, and is therefore 
taken as an estimate of the error introduced by the merging. At each step the 
vortices are examined pair by pair and the merging done only if the estimate 
is within a tolerance Vo- The exact estimate used is: 

|r i r 2I |*i - z£ 


<v Q 


(42) 


|r 1 F2I (Do + di) 1,5 (Do + d 2) 1,5 

where d\ and are the distances from z\ and Z 2 to the nearest wall and Dq 
is a parameter. The expression in Eq. 42 has the dimension of a velocity 
and is our estimate of the disturbance that a merging would impose on the 
boundary layer. Typically, Vo is of the order 10~ 4 Uoo or less. 

The disturbance estimate is the product of two factors. The first factor 
depends only on the circulations: 

. lEiEaL (43) 

|ri + r 2 | v 

Clearly, merging of vortices with large circulations is discouraged, as is the 
merging of two vortices that have nearly opposite circulations. (In that case 
the new vortex would be very far from the original ones. See Fig. 10b) 

The second factor depends only on the positions: 

|*i ~ *2l 2 


(Do + d I ) l - 8 (A> + «fe ) 1 - 5 


( 44 ) 
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Clearly vortices are more likely to merge if they are close to each other and 
far from the body. The parameter Do controls the relative variation of the 
estimate as a function of d\ and d, 2 - If Do is large (Do -\-d) has a slow relative 
variation near the wall and the density of vortices will be quite uniform. If Do 
is sn? "II (Do + d) gets very small near the wall, which discourages merging 
and will result in more small vortices subsisting near the wall. Thus, the 
parameter Do allows the user to shift the resolution from the wall region to 
the wake or vice versa, as illustrated in Fig. 11. 

The value of the tolerance Vq is not held constant for every time step. 
Instead, the programs slowly adjusts it to keep the number of vortices near 
the chosen number, raising Vo to make mergings easier if it sees too many 
vortices and lowering Vo if it sees too few. The number of vortices thus 
remains close to the input value, which is very desirable from a practical 
point of view. 

In contrast with some earlier approaches [22], this method of merging vor- 
tices is totally automatic and has a mathematical rather than a physical basis. 
For instance no effort is made to "manually" obtain a well defined Karman 
street; the vortices will probably take on this pattern at some distance from 
the body, but it will be destroyed as they move farther downstream. 

e) Numerical diffusion. 

The Vortex Method, at least in an unbounded fluid, has been shown to 
converge to the solution of the Euler (inviscid) equations. In reference [45] , 
(written with Dr. Y. Nakamura and Dr. A. Leonard) we applied the Vortex 
Method to several simple problems and by comparison with the known exact 
solutions confirmed the mathematical estimates: second order convergence, 
in terms of the core radius, was observed. These flows were all inviscid and 
unbounded, and the initial data had to be sufficiently smooth. Gaussian cores 
were used, but any core for which 7 is smooth and everywhere positive should 
also give second order convergence [19]. 

Although the Vortex Method solves the inviscid equations, there is evidence 
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of a significant "numerical”, or "parasitic”, diffusion in tlie solutions it 
produces. Essentially, strong velocity gradients induce strong accelerations 
which deteriorate the accuracy of integration of the motion of the vortices. 
This effect is different from the numerical diffusion present in many finite- 
difference methods, in that it is caused by velocity gradients instead of the 
velocity itself. Naturally this diffusion tends to zero as more vortices are used 
and shorter time steps taken, but it cannot be ignored when doing computa- 
tions with a practical level of resolution. To describe this numerical diffusion 
better we shall consider some of the invariants of the system. In this section 
only unbounded flows will be considered; the presence of a solid and the 
creation of vorticity at its surface would only complicate the discussion. 

It is well known that a system of point vortices is a Hamiltonian system 
[47] . The Hamiltonian of the system is the Kirchhoff function, defined by 

w = ]£ -sjr lQ g(ki - *A) (45) 

and the equations of motion become 


dxj _ dW 
* it - dy { 


(46.1) 




(46.2) 


Naturally W itself is an invariant. Other invariants are the total circulation, 
the first moment of the vorticity (equivalent to the momentum of the fluid) 
and the second moment (equivalent to the angular momentum) [14], [48]. W 
is also the energy of interaction of the vortices (their internal energy is infinite 
and has been separated from the interaction energy). These quantities are 
also invariants of the exact inviscid solution; these built-in invariants provide 
a basic advantage in using the Vortex Method. 

For a system of vortex blobs a Kirchhoff function can still be defined, by 
replacing the logarithm in Eq. 45 by an appropriate regular function. For 
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example Core 2 results in the Kirchhoff function 

lv =E^ i| °s(\/i r '- r ii 2 +‘ 72 ) 
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( 47 ) 


and Eq. 46 is satisfied. 

If we give all the vortices the same core shape, the vortex blob system has 
the same 4 invariants as the point vortex system (this is the ” semi-discrete” 
system: discretized in space but continuous in time). 

On the other hand it is clear that errors in the integration of the ordinary 
differential equations, Eq. 7, tend to scatter the vortices and therefore act in 
the same sense as a diffusion term. One way to describe this diffusion more 
precisely is to determine which of the invariants we mentioned are actually 
conserved and which ones are not in the solution of the fully discrete system 
(discretized both in space and time). Any deviation will be a result of the 
time integration errors. 

The total circulation is obviously conserved because each value T, is kept 
constant. The first moment of vorticity, / / ur, is conserved too if the time 
integration scheme is linear (which is the case for all the classical schemes) 
because it is a linear combination of the u's. 

The second moment of vorticity, / / cur 2 , is not conserved in general. It 
reflects the scattering of the vortices. In Ref. [45] we defined an effective 
viscosity u e by: 


4-00 4-o° 

/ / u(x, y)dzd$n 

— oo — oo 

(48) 

It is shown in Reference [48] that an exact viscous solution satisfies Eq. 48 
with u e replaced by v\ this motivates the definition of u e . The viscous 
diffusion steadily increases the second moment of vorticity, which is a measure 
of the spreading of the vorticity. The effective visccoity was calculated in 
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several cases and shown to tend to zero with the size of time step, for a 
given space discretization. This why we called the Vortex Method ”semi- 
inviscid” , meaning that the space discretization itself was not responsible for 
the diffusion, but that it allowed the time integration scheme to introduce a 
diffusive error (the same terminology is used when a method is said to ” semi- 
conserve” energy, i.e. it would conserve energy if the time integration were 
exact). 

The concept of effective viscosity according to Eq. 48 however has several 
weaknesses. It breaks down if the total circulation is zero, which is often 
the case, and if there are solid walls boundary terms appear which cannot 
be defined very reliably in a vortex simulation. More importantly, it is a 
global quantity. A concept that would yield a local effective viscosity would 
be much more useful, but has not been found yet. Thus it is not possible 
to produce the "modified equation” the way it is commonly done with finite 
difference methods, or to produce an "effective Reynolds number” of the 
computation. If that were possible, one could think of using the integration 
errors to ’ntroduce the desired diffusion. 

The reason why the second moment is not conserved is that it is not a 
linear combination of the r,’s; similarly, W is non-linear and will not be 
conserved. Delcourt and Brown used W (interpreted as the energy) for their 
definition of an effective viscosity [31]. The effective viscosity turned out to 
be positive, since the energy decayed steadily. For the time integration they 
used the Euler explicit and the Huen scheme (also called Runge-Kutta, first 
and second order). 

We showed in Ref. [40] that the non-linearity of Eqs. 7 and has a strong 
influence, even in a very simple case: two vortices isolated >n space. If their 
circulations are Ti and r 2> their (complex) positions are given Z\ and Z 2 , and 
they are treated as point vortices, then their spatial separation Z — (Zi — Z 2 ) 
satisfies the first order ordinary differential equation: 

dz .(r,+r 2 ) 



If Zo is the initial separation, at t = 0, the solution is: 

Z(f) = Z 0 e’ n °‘ (50) 

where Oo is defined by Oo — (ri + r2)/(27i|Zo| 2 ). The vortices orbit 
together, and the second moment of vorticity is constant; in that sense the 
discretization by vortices did not introduce any diffusion. The linear ordinary 
differential equation: 

^ = in 0 Z (51) 

has the same solution and is more familiar. In finite difference methods, 
convection terms generally produce pure imaginary eigenvalues, which makes 
Eq. 51 a good model problem. 

Although the two equations have the same exact solution their numerical 
integration, by the same scheme, can give widely different results in terms 
of stability. We shall concentrate on the modulus of Z since we are mostly 
interested in scattering. Fig. 12 shews |Z| as a i unction of time as found 
in a numerical solutions to Eqs. 49 and 51 for a typical case: Adams- 
Bashforth 2 and Lomax schemes, and several values of the time step. The 
Lomax scheme is especially adapted to the integration of Eq. 51 [46] . With 
Adams-Bashforth-2 the nonlinearity of Eq. 49 reduces the error, compared to 
the linear equation, because as |Z| increases the angular velocity decreases, 
and the integration becomes more accurate. The solution to Eq. 51 grows 
exponentially, which is a strong instability, while the solution to Eq. 49 only 
grows like f 1 / 8 . The integration of the linear equation by the Lomax scheme 
shows exponential decay, while the integration of the non-linear equation 
keeps |Z| finite. In this case the second moment increases and decreases 
periodically; the effective viscosity is not constant, and even takes on negative 
values. It appears that the angular velocity cannot remain below a given value 
(about 0.24/Af with the Lomax scheme). 

Similarly, in practical cases with many vortices the time integration scheme 
does not allow values of angular velocities above a certain level, and scatters 
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the vortices when the local velocity gradients are too stiong (in that sense 
the Vortex Method is stabilized by its own integration errors). 

In the computations done in this study, a typical value of the angular 
velocity in the near wake is 10, and the time step is 0.03 (in compatible 
units). The product is non-dimensional and of the order of 0.3: obviously 
integration errors will be significar*. at the scale of the individual vortices. 
Considering the fact that at high Reynolds numbers the angular velocities in 
the exact solution would be as high as several hundred units, at least an order 
of magnitude greater, it is also clear that resolving all scales is not possible. 
The Vortex Method performs well in spite of such errors partly because the 
conservation of circulation and momentum are built-in. 

The effect of merging can be examined in the same spirit. When two 
vortices of the same sign (the more likely case) merge the second moment of 
the vorticity distribution decreases. The merging concentrates vorticity and 
this is especially apparent in the far wake. If the mergings occur far enough 
from the solid body the effect of this "reverse diffusion" is small. 

2) Inner flow. 

The boundary layer equations, Eqs. 29, 30 and 31, are solved using a 
finite difference discretization in space and an implicit method in time. The 
accuracy is of second order in one space direction, fourth order in the other, 
and first order in time. The Baldwin-Lomax algebraic turbulence model is 
used. 

a) Extent of region and discretization in space. 

The region is a band of width 6 placed around the body. 6 is small, com- 
pared to the radius of curvature, and the curvature of the band is neglected. 

The functions u, u and v are assigned values at the nodes of a grid. The 
grid is stretched in the s direction, according to the distribution of the points 
along the wall. In the n direction an exponential stretching is used to give 
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a finer resolution near the wall. This is especially useful for turbulent cases, 
in which the viscous sublayer is very thin. With the grids used the cell size 
near the wall was of the order of 5, in "wall units”, so that the first point 
was well within the viscous sublayer [50], The equations are transformed to a 
computational plane where the grid is uniform. Centered differences are used 
for all derivatives in all directions. Second order accurate differences are used 
in the n direction, in which the grid can be made very fine without penalty. 
In the s direction, the grid repoduces the intervals involved in the outer flow, 
and these cannot be made very short. For this reason, and to make the 
convection of vorticity as accurate as possible, fourth order accurate Pade 
differences are used in the s direction. Naturally, to actually obtain fourth 
order accuracy the grid should be smooth enough, which is not always e*. y 
when generating complex shapes. 

b) Computation of the velocity field. 

The u velocity in the grid is obtained by integrating Eq. 30 upwards from 
the wall, where u = 0. The v velocity is then obtained by integrating Eq. 29 
with v = 0 at the wall. In both cases, the second order accurate "trapezoidal 
rule" is used in the n direction, and du/d$ is obtained by Pade differences. 

c) Time integration scheme. 

The integration in time is done using a first order accurate implicit scheme: 
it is the Euler implicit scheme, except that the velocity components are 
"frozen” at the old time level. Completely linearizing the non-linear convec- 
tion term, U.Vcj, would make the matrix inversion much more costly without 
formally improving the accuracy. Furthermore the Euler implicit scheme is 
very stable according to a "frozen velocity" analysis, and there is no evidence 
that the incomplete linearization hurts its stability. 

This first order scheme is used because the time variations in the boundary 
layer are very slow (typically 300 steps per period) and because implementing 
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a higher order scheme would be much more complex, again because of the 
difficulty in linearizing the velocity components. 

The boundary condition in the s direction is periodic and does not require 
special attention. The boundary conditions in the n direction are implicit, 
which is desirable in the presence of a viscous term and with a fine grid. 

d) Approximate factorization. 

The time evolution equation is written in "delta” form and the implicit 
operator is approximately factored into two tridiagonal operators, one in each 
direction. This simplifies the solution without degrading either the first order 
accuracy in time, or the accuracy of the steady state. 

The operator in the s direction is periodic and tridiagonal. It is solved 
by the Thomas Algorithm, adapted to periodic matrices. The operator in 
the n direction has the three diagonals plus a full first line representing 
the integral across the layer. This integral condition replaces the boundary 
condition at the wall; more details will be given in the chapter on the coupling. 
The boundary condition at 6 will also be described later; it is included by 
modifying the last line of the matrix. This matrix is solved by the Thomas 
Algorithm, this time adapted to start the elimination from the bottom and 
eliminate the first line too. 

e) Artificial dissipation. 

Finally, an artificial dissipation is added in the s direction. The centered 
differences used to approximate the first derivatives in the s direction do 
not couple the even and odd lines, and a small positive term representing a 
derivative of even order is added to the time derivative to absorb energy and 
avoid the appearance of oscillations. A fourth order derivative is generally 
used, to disturb the slow varying components as little as possible. Depending 
on the amount added, it might be necessary to treat the artificial dissipation 
implicitely to preserve stability. 
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Probably because of the constraints imposed by the outer flo-w, and of the 
non-linearity, a catastrophic instability does not result if the dissipation is 
omitted in this program: the solution remains bounded. However it exhibits 
short wave undulations in the s direction, which have half-periods close 
to the grid size, and are very probably caused by the inaccuracies of the 
finite difference method when treating the convection term with in '-constant 
velocity. Naturally the coefficient of artificial dissipation, t, is given a value 
as low as possible. This will be illustrated in the "Results” chapter. 

As a whole the numerical method used for the inner flow closely follows the 
theories developed by Beam, Warming, Steger and Pulliam at NASA Ames 
[49]. 

f) Turbulence model. 

The Baldwin-Lomax turbulence model was chosen because it is simple to 
use and was designed for separating flows [50] . It is based on the Cebeci- 
Smith model, but modified to ensure a normal behavior even when the 
boundary layer thicknesses become very large. It is an algebraic, or "zero- 
equation” , model; it does not require the solution of any additional differential 
equations, or any special conditions at the outer edge of the inner region. 

The turbulent stresses predicted by the Baldwin-Lomax model are multi- 
plied by an intermittency factor £ which is a function of s only, and switches 
from 0 to 1 as the boundary layer undergoes transition. The transition 
model proposed by Baldwin and Lomax is not used; it does not seem to 
take sufficiently into account the pressure gradients which are very strong in 
flows around cylinders for instance. The criterion described by Schlichting, 
which is based on his own stability theory and an empirical correlation by 
Granville, incorporates the dependence on the pressure gradient and is used 
instead [51] . This transition model produces a position s t of instability of 
the laminar boundary layer, and a position of full transition. This delay 
is used in the program to make the transition smoother: the intermittency 
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factor fi(s) is defined by: 


fi{s) = 0 if s < 5, (52.1) 

P(s) = 3^ 3 ~ if 5, v s < s t (52.2) 

P(s) =1 if 5 > St (52.3) 

The choice of a cubic function for Eq. 52.2 was arbitrary; it was chosen 
merely to make ft a smooth function of s. 

The turbulence model is present in the algorithm regardless of the Reynolds 
number. However, for Reynolds numbers of 10 5 or less, transition is not 
predicted (although instability often is) and the turbulent stresses are never 
activated. Thus the computation is fully laminar in such cases. 

3) Coupling. 

The interaction of the two regions involves the matching of the velocity 
fields and the transfer of vorticity acrr , the interface. 

The velocii.es are matched by properly setting up the vortex sheet that 
represents the inner flow vorticity in the Biot-Savart integral. The outer 
velocity field is a function of the r/s and IYs, the values of U, and the position 
of the vortex sheet. If it satisfies Eq. 10.4 then the velocity at the wall (under 
the sheet) will be zero, and the tangential velocity over the vortex sheet will 
be U,. Thus U e is the value of the tangential outer velocity at n = 6. 

On the other hand the tangential component of the inner velocity at n = <5 

is: 

6 



n— 0 


in view of Eq. 30. For the tangential velocities to match, this integral must 
be equal to U e . 
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Matching the normal velocities is not as crucial, because the normal velocity 
is still small at 6. The program runs quite well without any effort to match 
the normal velocities. However some finer effects can be added by doing so. 
Since the vortex sheets is a simplified representation of a layer of vorticity of 
finite thickness 6, the vertical position of the vortex sheet is arbitrary within 
the thickneu 5; it is natural to place it at the centroid of the inner flow 
vorticity, defined by: 



Jn, o nu(s,n)dn 
fn= 0 u{s,n)dn 


(54) 


Fig. 13 shows that if the vortex sheet is placed at 5* the normal velocity at 
a distance 6 from the wall must be 




(55) 


for mass to be conserved. 

Now from Eq. 29, the normal component of the inner velocity, at n = 6, 
is: 


v 

”(M) = - J 


du 

ds 


dn 


(56) 


n— 0 


Using integration by parts we obtain: 


v(M) = - 


du 

i"s7 


10 


Jo 


~ Js I n ^ dn 


(57) 


n~(J 


This equation, combined with Eqs. 53, 54, and 55, and the fact that 66 /ds — 
0 shows that the normal velocities at 6 indeed match. 

The vortex sheet is placed at 6* because 5* is the centroid of the vorticity. 
However integration by parts shows that, if the boundary layer is entirely 
contained in the band of thickness 6, 6* is the classical displacement thick- 
ness defined in boundary layer theory [51] . The line defined by <5* acts as an 
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effective boundary for the outer flow; its slope introduces the small amount 
of normal velocity due to the thickening of the boundary layer. This yields a 
boundary layer procedure of higher order and the incorporation of the dis- 
placement effect is necessary to avoid the singvdarity that otherwise appears 
in the boundary layer solutions near separation or reattachment [52] , [53] , 
[54l The boundary layer acts on the pressure both by the emission of vor- 
ticity into the outer region and by the displacement effect inherent in 5*; the 
outer region vorticity controls the broad features of the pressure distribution 
while the 6* effect corrects the pressure locally, especially near separation. 

Although many questions are still open the consensus seems to be that a 
truly unsteady boundary layer, even in a direct solution ( U t imposed), does 
not have a singularity at separation [52] , [53] , [54] , [55] , [561 . However, if 
the solution as time progresses tends to a steady state, the shear stress dis- 
tribution will steepen and tend to a singular d>st r t button unless the boundary 
layer is allowed to relieve the pressure gradient by the displacement effect. 
This is what is sought in this algorithm; fairly smooth solutions are obtained 
but some oscillations near separation suggest that the problem might not be 
entirely solved. Naturally, the production of fair numerical solutions is not 
a proof of the regularity of the differential system unless a thorough conver- 
gence study is made like in reference [54] . This was not possible here, mostly 
because of the high cost of the vortex part of the computation. 

In some plots of the computated results (especially in Figure 36), the vortex 
sheet is shown as a solid line over the wall and it is apparent how it lies 
very close to the wall in the part of the boundary laye* having a favorable 
pressure gradient, then leaves the vicinity of the wall, until the boundary 
layer separates and injects itself into the outer region, becoming a free shear 
layer made of vortices. 

One problem persists regarding the positioning of the vortex sheet at <5*. It 
is that 5* can take on negative values, or values larger than 5, and in general 
is not very smooth in the regions where U e is small. Values of <5* larger than 
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8 are not acceptable, because then some vortices would be under the vortex 
sheet and their tangential velocity would not be correct. In addition, we hare 
seen that the vortex sheet should be kept very smooth. For these reasons, 
the function 8* is filtered and truncated between 0 and 8, giving 8, and the 
sheet is placed at a distance 8 over the wall, instead of 8 . 

The thickness 8* is obtained from the inner flow solution. The determina- 
tion of the quantity U e is more complex and coupled to the transfer of vor- 
ticity. The whole procedure will now be described; it reflects the flux of 
information from one module to the other aud was devised on an intuitive 
basis. The flow chart in Fig. 14 illustrates it. 

A buffer is used that communicates alternatively with the outer flow and 
with the inner flow. It is a vortex sheet of intensity Bf. Starting from the 
outer flow, at each time step the vortices that crossed the interface are put 
into the buffer and considered as candidates for absorbtion by the inner flow. 
The buffer then communicates with the inner flow. 

The circulation per unit length, under the old vortices, will be (U e -f Bf). 
Equations 10.3, 10.4 and 10.5 are then solved. In one time step the vorticity 
that is generated 'at the wall will not reach the outer flow in significant 
amounts (the vorticity diffuses to distances of order y/iTKi and vVA t will be 
.002 or less while 8 will be .015 or more). Therefore Eq. 10.3/10 4/10.5 can be 
solved to a very good approximation by considering the outer flow vorticity 
as known and the strength (U e + Bf) of the inner shear layer as unknown. 
This amounts to assuming that for one time step the flux of vorticity through 
the wall, which is also the pressure gradient, does not depend on the shift of 
inner .egion vorticity in the n direction. 

In that sense the boundary layer solution is "direct” at each step and the 
coupling is not "strong” in the sense of [53] ; the pressure distribution is 
imposed on the boundary layer for this step, and will respond to the boundary 
layer only for the next time step. This should be sufficient since the variations 
in the boundary layer are very slow. 
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[U e -f Bj) is computed by solving a linear system. If N w is the number 
of wall intervals, the N w unknowns are the values of (U e + Bj) over each 
vs 11 interval. The first N w — 1 equations govern the differences between the 
V3 ! ues of the stream function at the N w wall points (this is the integral form 
of Eq. 10.4, which is considered less sensitive and therefore more efficient 
th in the colocation form), and the last equation governs the total circulation 
ea itted by the solid (Eq. 10.5). The stream function is the sum of the 
sti earn function generated by the freestream and the existing vortices (which 
is tnown and forms the right hand side), and of the stream function generated 
by the vortex sheets (which is the unknown). The matrix is computed at the 
beginning of the run and "Gauss eliminated” once and for all, since it does 
not change. Note that a rotation of the solid does not affect the matrix; but 
if several solids were to move independently, this would affect the distance 
between wall points, and a different matrix would have to be inverted at each 
tine step. The cost would then be higher. 

The linear equations are set up to strongly couple each unknown with the 
eqiation of the same index, so that the matrix has its larger elements near 
thi diagonal. As a result, the matrix is well conditioned enough for Gaussian 
elimination with partial pivoting or even without pivoting (both on the CDC 
7610 and the CRAY). 

Then the inner flow is advanced. In particular, the transfer of vorticity 
th 'ough the interface is computed. This vorticity is transferred between U e 
and B / , but (U t -\-Bf) takes the smlue that was just computed. Two cases are 
possible: v e < 0, and v e > 0. v e is the normal velocity at 8 and is obtained 
frcm the inner solution (Eq. 56). As was done for the inner flow solution, the 
ve ocity v e lags by one time step. 

r i v e < 0 the transfer is imposed by the outer region, in keeping with Eq. 
30 The buffer vorticity is injected into the boundary layer; after the transfer 
tli ; buffer is empty and all the vorticity is in U e . This injection of the buffer 
coistitutes the boundary condition at n — 8 for the inner flow. 
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If v e > 0 the transfer is imposed by the inner flow. The inner flow rejects 
the buffer and injects additional vorticity into it. The flux is v e u e , where w e 
is the value of the vorticity at n — 6. This in itself constitutes the boundary 
condition at 6, since it amounts to setting the viscous flux to zero (it can also 
be interpreted as a linear extrapolation). 

In both cases, the sum of Bj and of the integral in Eq. 53 is equal to the 
value that was computed for ( U e 4- Bj). This yields the integral condition 
for the inner flow vorticity. 

The new values of U e and Bj have now been determined and the program 
returns to the outer region. The U e vorticity stays in the boundary layer and 
U e gives the strength of the vortex sheet. The buffer vorticity is injected into 
the outer flow under the form of new vortices if v e > 0. (If v e < 0 the buffer 
is empty.) The values of 5* have also been computed and the vortex sheet is 
repositioned. 

The outer flow is then advanced, which involves the computatfon of the 
velocities, the motion of the vortices, and the mergings. The program is now 
ready to start a new loop by determining the vortices to be put into the buffer 
(Fig. 14). 



IV) RESULTS. 


1) Parameters. 

The main parameters are the number of vortices, N v , the time step, At, the 
core radius, a, the distance, Rq, from the wall where newly created vortices 
are placed, the grid thickness, 8, and the artificial dissipation coefficient, e, 
if KPD3 i, c used, and the parameter, Do, of the merging device. 

a) Number of vortices, N v . 

The cost of a computation depends strongly on N v since the computer time 
per time step is roughly proportional to N„; naturally, the larger N v is the 
greater the details that are reproduced and the more accurate the simulation. 
It is impossible to rationally select a minimum value of N v for a particular 
situation. One should observe the solution and ascertain whether the smallest 
features considered significant contain at least several vortices and therefore 
have some degree of structure and some ability to be strained. (A good 
graphics system is essential for the monitoring of vortex computations). If 
an eddy contains only one computational vortex ir is effectively circular, and 
of size o. 

A more quantitative estimate of accuracy, in selecting the value of N v to 
be used, is the tolerance V 0 for merging. The larger N v is the later the 
mergings will occur and the smaller V 0 will be. As an exam ole the same code 
(XPDl) was run for the same case (a square) with N v = 800, then with 
N v = 1200 (F ; g. 15). At the end of the simulations Vq was 7.6 x 10~ 4 in 
the first case and 1.7 x 10“ 4 in the other case. The difference is significant: 
with 800 vortices, mergings occured that caused a disturbance four times as 
strong as would be allowed with 1200 vortices. However both values are small 
compared to Uoq. 
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For a single smooth body a value of 1000 is generally sufficient: the 
significant eddies in the flow are not very small and the other sources of error 
(boundary layer assumptions, turbulence model, etc.) probably dominate. If 
the body has sharp corners or a trailing edge then larger values of N v might 
be desirable. The need for resolution is also stronger if the separating boun- 
dary layer is very thin; for example, the circular cylinder in the critical range 
of Reynolds numbers. In that case 1600 vortices were used. Finally, with 
several bodies the wake of the first body interacts with the other bodies; then 
it is justified to use much more vortices. The Vortex Flowmeter simulations 
described in Reference [41] used N v = 3600. 

b) Time step, Af. 

As was the case for N v , the choice of Af is a compromise between cost and 
accuracy. 

The Lagrangian method can be quite accurate in the wake without a very 
short A t because accelerations are moderate there, which makes Eq. 7 easy 
to integrate accurately. Similarly the Eulerian method in the boundary layer 
can be accurate without a very short step: the boundary layer often evolves 
slowly and in that case C. F. L. numbers much larger than 1 are acceptable 
(the C. F. L. number is |[7|At/As, with U the local velocity and As the grid 
size in the s direction). The region that demands a short step is in general 
the intermediate region, and this is for two reasons. First, the vortices that 
are just outside the inner region often pass several grid points in one step 
while the inner region points interact only with their immediate neighbours. 
This can create an imbalance because the signals do not travel at the same 
celerity in the two layers. Second, if the time step is long the newly-created 
vortices are stronger (their circulation is proportional to Af) and such strong 
vortices disturb the inner region. 

To estimate an acceptable value for Af the user should observe the simula- 
tion where the body has tight curves; if Af is too long the vortices will not 
follow the wall. The C. F. L. number should not be much larger than 1. 
In general the results depend on At more directly than on N v \ one of the 
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weaknesses of the method is that it is only first order accurate in terms of 
At. 

c) Core Radius, a. 

Unlike N v and At the core radius does not influence computing cost and an 
optimum value exists, instead of a compromise between cost and accuracy (in 
Ref. [45] we systematically determined this optimum value in a few simple 
cases). If cr is large, the velocity is very smooth locally and the noise is low: 
as a result the vortices will not scatter much. On the other hand a large 
core radius can suppress velocity gradients that are physically significant and 
"freeze” a coherent structure that would be better represented if the cores 
were small enough and allowed it to evolve. 

Fig. 16 shows the same flow computed with a — .005 and then a = .05. It 
is clear that the simulation is not very sensitive to the value of <r. changing it 
by a factor of 10 did not cause a striking difference. As a rule, a should be of 
the order of As/2, where As is the spacing of the points along the wall. The 
value of <7 influences the coupling between wall points and creation points, in 
the same way as the value of Rq does (see subsection d)). 

d) Distance Rq. 

The points where new vortex blobs are introduced are located at a distance 
Rq over the wall. In addition, vortices that are found within a distance Rq of 
the wall are treated as being absorbed by the wall layer. Thus Rq is a rather 
important parameter. 

If KPD3 is used, Rq is equal to 6, so that the vortices are created at the 
edge of the viscous region. If Core 1 is used, Rq is equal to the core radius 
< 7 , so that the edge of the core is tangent to the wall. 

If Core 2 is used in KPDl or KPD2, Rq is an independent and non- 
physical parameter. A good value for Rq is aoout As/2, where As is the 
spacing of the points along the wall. Much smaller values would let the 
vortices go too close to the wall points (where the stream function is sampled) 
and create noise in the pressure. Much larger values would weaken the 
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coupling between each creation point and the wall point below it; again, the 
result might be oscillations in the pressure distribution. Such oscillations 
are a sign that the system is not functioning properly, and will also strongly 
disturb the integral boundary layer solver. 

For simplicity Rq is held the fixed for all the points along the wall. On 
the other hand As might vary, for instance if wall points are clustered in a 
region that is thought to require more resolution. (The selection of the wall 
points and their clustering is left to the user). To keep the ratio Rq/As at a 
value of the order of 1/2, the clustering of points should be moderate. In the 
applications presented here points were clustered near sharp edges or trailing 
edges, or sometimes in the separation regions, but the ratio of the largest 
value of As to its smallest value did not exceed about 2. 

e) Selection of 6, for KPDZ. 

The thickness 6 of the computed viscous region results from a compromise 
and can be chosen by observing the solution. 

On the one hand, the larger the value of 6, the greater the domain treated 
by the viscous solver, which is good. (In addition, extending the computa- 
tional viscous region is not very costly). The attached part of the boundary 
layer clearly must be contained in the grid; the vorticity contour plots are 
helpful in ascertaining this. Another i?ay to assure it is to compute the dis- 
placement thickness, 6*: it should be of the order of <5/2 or less along the 
entire attached boundary layer. Near separation it is normal for 5* to be- 
come comparable to 6 or even exceed it (vorticity of the two signs is present, 
and as a result the centroid can have large excursions). Fig. 17a shows a 
computation in which 6 was too small: the boundary layer reaches the edge 
of the grid, even in a region with favorable pressure gradient. This defeats 
the purpose of having an inner viscous region. 

On the other hand, the larger 6 is the less the boundary layer assumptions 
are justified. The errors associated to these assumptions grow. These errors 
are hard to estimate quantitatively, but the results often give indications 
when 6 is too large: the stability decreases and oscillations appear near 
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separation. Probably, the variations of the displacement thickness become 
too steep and disturb the algorithm. Fig. 17b shows such a case. 


f) Selection of e. 

A small amount of artificial dissipation proves to be necessary to keep the 
solution smooth in the inner region, so that the finite difference approxima- 
tions are accurate. Since fourth-order dissipation is used, it is difficult to 
compare the effect of the artificial dissipation with other sources of dissipa- 
tion, for example the viscous stresses which are a icond-order term. Fig. 
18a shows a simulation with e = 0. Ripples appear in the vorticity contours 
and the other quantities involved with the inner region. The ripples are in the 
s direction, which was to be expected since that is the direction in which the 
grid is coarse and viscous stresses comparable in magnitude to the convection 
terms are not present. 

This is why artificial dissipation is added only in the s direction. In Fig. 
18b the flow is simulated with the grid and other parameters the same, but 
e =ss 0.8 x 10“ 5 . The solution is now smooth, e can be increased further, 
even by a factor 10, without any apparent effect. This is important since it 
shows that the dependence of the solution on e is weak, provided that e is 
large enough to eliminate the ripples. 


g) Merging parameter, Dq. 

The effect of Do was described in Chapter HI. In this case also, the only 
way to determine a good value for the parameter is to observe the simulation. 
However, a good rule is to make Do about 5% of the body size if only one body 
is present. If there are several bodies, more vortices are devoted to computing 
the wake of the first body, so that it is not too coarsely represented when it 
strikes the second one. To achieve this, D 0 is made larger: about 50% of the 
distance between bodies. 
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2) Memory requirements and computing times. 


a) Codes KPD 1 and KPD2 

The codes KPD 1 and KPD2 have very similar requirements; the integral 
boundary layer solver requires very little memory and computing time. The 
total memory is about 10 5 words, and little effort was spent trying to reduce 
it. It could certainly be reduced to about 0.7 x 10 5 . With 200 wall points the 
matrix alone requires 0.4 x 10 5 words. The computing time with N v = 1100 
and N w = 200 is about 0.4 seconds per step on the CRAY-1. This figure is 
probably close to the minimum; the program was carefully written to reduce 
the CPU time and all the major operations were vectorized. Theses operations 
are the computations of the interactions (O(iV^)), the merging tests (O(N^)), 
and the computation of the stream function at the wall points {0(N v N w )). 
The non-dimensional time step UooAt/c was 0.015. Thus a complete period 
of oscillation for the square body requires about 250 seconds on the CRAY-1. 


b) KPD3 code. 

The KPD3 code requires a fair amount of additional memory for the 
differential boundary layer solver; the total is now about 2.8 x 10 5 words in 
the high resolution runs (again, no special effort was devoted to lowering the 
memory requirements). The additional CPU time required for the boundary 
layer is modest: about 15%. The other operation that is more time-consuming 
is the computation c' » '.o velocity field induced by the vortex sheet at the wall. 
It is 0(N v N w ) and now involves a complex logarithm, whereas in KPDl and 
KPD2 only arithmetic operations were involved. 

Simulations of the circular cylinder with KPD3 were run with two levels 
of resolution. At the lower level N v = 1100, N w — 200 and At = 0.025, and 
the CPU time was about 1.1 seconds per step, or 450 seconds per shedding 
cycle. The high resolution runs used N v = 1600, N w = 300 and At = 0.02, 
resulting in a CPU time of 2.1 seconds per step, or 1050 seconds per cycle. 


68 



3) Results from the program KPD 1. 


a) Flow around a square. 

The flow around a square is a good test because the viscous effects seem to 
be reduced to a minimum. The experiments exhibit a flat drag curve from 
Re = 10 4 to Re = 10 7 [57], probably because the primary separation occurs 
at the front corners irrespective of boundary layer thickness or of its laminar 
or turbulent character. 

The computed mean drag coefficient is 1.8; the experimental value is 2. 
The computed Strouhal number is .11, which agrees with experiments. These 
results were obtained with Core 1, N v = 1000, N w = 320, a = Rq — .02 
and At = .03 (C/oo = 1 and c = 2). The tendency to underpredict the 
drag of the square is real: with other sets of parameters the drag was often 
even lower than 1.8, of the order of 1.4. This is a little disturbing since the 
case of a square was chosen precisely because complex viscous effects are not 
thought to be present. One possible cause for the inaccuracy is the difficulty 
for the algorithm to model the flow near a sharp corner, where the radius 
of curvature is small even compared to the size of the vortex cores and the 
spacing of the vortices. 

b) Starting vortex at a sharp corner. 

This example illustrates "viscous” behavior in the simulations done using 
the Vortex Method. In this case the viscous character is properly, although 
fortuitously, reproduced. 

The body is a diamond and is a fair representation of the experimental 
situation in Ref. [1] , for short times after the start. The case chosen here 
is an angle of 6£P and a constant velocity after the impulsive start = 1/3 
and m = 0 with the notation in [58]). We shall focus on the flow pattern 
near one of the corners a short time after an impulsive start (the Kaden 
Vortex), and compare it with the experimental pattern. At time zero the 
flow is irrotational and as a result the velocity tends to infinity at the corner. 
This is the solution to the inviscid equations and it would stay unchanged 
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at all times. However, in the numerical solution, vortices are seen to leave 
the wall and form a starting vortex that resembles a spiral (Fig. 19a). The 
streamlines also show that the flow leaves the surface at the corner along the 
tangent to the upstream face; this pattern is generally accepted as correct 
and consistent with the Kutta condition. The streamlines are also shown 
in Fig. 19a and reveal the topological structure of the velocity field, with a 
half-saddle point at, the corner and one on the back face. 

Visual agreement with experiments is good. In addition, the growth of the 
length scale with time was compared with similarity theory. This theory, for 
an angle of 60°, predicts that the length scale will be proportional to t 5 / 7 
(B. Cantwell, personal communication). In Fig. 19b, the distance from the 
tip to the point of zero velocity (near the center of the spiral) is plotted as 
a function of time, in Log-Log coordinates. The curve is close to being a 
straight line and its slope close to 5/7. The Vortex Method has succeded in 
simulating the formation of a starting vortex and associated establishment 
of the Kutta condition. 

In this case it is clear that the motion of the vortice , being advanced 
finite step by finite step, cannot accurately follow the sharp kink in the wall, 
especially if the magnitude of the velocity is large. After the first vortices 
have separated from the corner the flow pattern changes, the streamlines 
become smooth and the velocities smaller, and the simulation can be quite 
accurate. However the process, originally, is caused by inaccuracies in the 
numerical solution of Eq. 7. 

c) Airfoil at ,ow incidence. 

This case is shown only to illustrate how KPDl, a pure Vortex algorithm, 
can become inaccurate when long stretches of attached flow are present. This 
is a major issue about the Vortex Method in general. 

The flow around an airfoil at a low incidence, 10°, is shown on Fig. 
treated by KPDl. The correct solution is a flow that renr.lr:> ;.A x be* . 
almost to the trailing edge, with a very narrow wake and ver low ui / I be 
lift is accurately specified by potential flow theory and the K<..tto ' ,ndi 
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is satisfied. 

In the numerical simulation, on the other hand, separation (defined as a 
departure of the flow from the solid surface) is seen to occur just downstream 
of the leading edge on the upper surface. As a reside the wake is much too 
wide, the drag is significant and the lift too low, probably because the trailing 
edge is immersed in the wake and the Kutta condition is not satisfied. 

The strong thickening of the boundary layer is a clear manifestation of 
the numerical diffusion that was described earlier (Chapter III). Basically, 
the boundary layer as represented by vortices can remain attached only in a 
significantly large and favorable pressure grad ent (it actually does along the 
lower surface). Thus KPD 1, and pure Vortex methods in general, are not 
well adapted to the simulation of flows past streamlined bodies. For these 
bodies, especially for airfoils, the accuracy will improve dramatically when 
the boundary layer is treated properly. In Fig. 20b the same flow is shown, 
treated by KPD2. The pattern is now correct and the quantitative features 
like lift, drag and moment much more accurate (see section 4). 

4) Results from the program KPD2. 

KPD2 is used mostly to compute airfoil flows. KPDl would not be 
accurate because airfoil flows have long boundary layers, and KPD3 is still 
restricted smooth shapes. 

a) Attached flow on an airfoil. 

This example is a direct test of the accuracy of the method, in an admit- 
tedly simple case, by comparing it to an exact solution. A Joukovsky airfoil is 
treated at an incidence of 5°. The solution for potential flow with the Kutta 
condition satisfied is known analytically; the pressure coefficient predicted 
by the numerical method together with the exact one are compared in Fig. 
21. The agreement is very good, which proves two facts. First, the flow as 
a whole is very well predicted and the circulation has been properly chosen 
by the algorithm to satisfy the Kutta condition. Second, the method of com- 
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putation of the pressure, using Eq. 13, is accurate; this is important because 
Eq. 13 does not involve Bernoulli’s theorem, which happens to be valid here 
but is not valid when vortical flows are simulated. 


b) Starting vortex on aa airfoil. 

Formation of the starting vortex and establishm ent of the Kutta condition 
are key features cf the flow around airfoils; they are responsible for the exis- 
tence of lift. Traditionally, the Kutta condition has been added to inviscid 
models in order to mimic a viscous phenomenon, namely the separation of 
the boundary layer that takes place at the trailing edge if the fluid attempts 
to flow around it. However the Vortex method, although it is in principle 
inviscid, reproduces this viscous feature without any intervention. This is 
another case in which the algorithm conveniently disobeys the inviscid equa- 
tions: vortices are "thrown away" from the trailing edge if large velocities are 
present, until the circulation is correct and the trailing flow is smooth. 

A convincing and probably accurate simulation of the starting vortex has 
been obtained and compared with the results of Wagner’s theory [2] . An 
airfoil is started impulsively in an irrotational fluid. The airfoil is at incidence 
but there is no circulation around it; thus the Kutta condition (in classical 
terms) is not satisfied at time zero. The boundary layer is seen to separate 
from the lower surface at the trailing edge in Fig. 22. The vortex sheet 
it carries curls up into a vortex that is swept downstream. Since the total 
circulation remains zero, a circulation of the opposite value is established 
around the airfoil. This circulation grows with time, as does the lift on 
the airfoil. This lift starts at about half the steady value; Wagner’s theory 
predicts an initial lift of exactly half the steady value. In general, agreement 
between the two curves, shown on Fig. 23, is very good. Wagner neglected 
the curling up of the vortex sheet and this might account for some of the 
disagreement. Also, Wagner used thin airfoil theory and the steady lift curve 
slope was 27r; here the airfoil was 12% thick, which results in more lift. The 
convergence of the lift to the steady value is made slower by the downwash 
of the starting vortex; this downwash decays only like t~ l . 
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Also of interest are the motion of the center of pressure, and the presence 
of drag, due to the downwash created by the vortex. At later times the center 
of pressure is at the quarter-chord point and the drag is zero (Fig. 22). 

c) Dynamic stall. 

The next case is the dynamic stall of the same NACA 0012 airfoil. Dynamic 
stall is a challenging problem, and a cate was chosen for which experimental 
results are available [60]. The airfoil performs prescribed oscillations in pitch 
with the pivot at the quarter chord. The incidence is a sinusoidal function of 
time given by 

a(f) = Qo -f- atisin(kUoot/c) (99) 

The Reynolds number is 2.510®. In the simulation, the lower boundary 
layer was in a favorable pressure gradient and generally it did not undergo 
transition and did not separate until the trailing edge. The upper boundary 
layer tended to undergo transition ana rer: sin attached at low angles of 
attack; at high angles of attack, it separated while still laminar. It switched 
instantaneously from one state to the other; this is not very satisfactory, but 
is inherent in the transition model that was used. Also, separation bubbles 
in which the separated shear layer makes a transition and reattaches on the 
upper surface cannot be reproduced by the algorithm. This i3 unfortunate 
since such a bubble is probably present, at least during pan of the cycle 
( KPDZ was written to eliminate these shortcomings as far as possible). 

The history of the flow durng one cycle can be followed on Fig. 24. The 
airfoil together with the vortices are shown; the force is displayed in terms 
of its magnitude, direction and axis of application (the center of pressure). 
Finally, the dashed line shows the suction distribution, measured normal to 
the surface and referenced to the stagnation pressure. The cycle begins with 
attached flow at 5° incidence; the ivutta condition is satisfied as evidenced by 
the smooth wake emanating from the trailing edge. The center of pressure is 
very close to the quarter-chord, which is expected for a symmetrical airfoil. 
As the incidence rises, counterclockwise circulation is shed; lift on the airfoil 
rises and the suction peak at the jose gets larger. At ?'«out 20°, tne pressure 
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gradient that follows this suction peak becomes too severe and the upper 
boundary layers fails to undergo transition; it separates instead from the nose 
of the airfoil. A vortical region of moderate thickness forms on the upper 
surface. The flow then reattaches for a short time before entering the fully 
stalled condition (it is not certain whether this reattaenment is physically 
correct). 

After the incidence passes its maximum of 25°, the flow on the upper surface 
remains separated and large eddies form on the upper surface. These eddies 
are accompanied by low pressure regions. The surface pressure distribution is 
strongly disturbed; the center of pressure moves away from the quarter-chord, 
mostly towards the trailing edge. The iift stays roughly constant until the 
hr ier s back to about 15°, then falls. During the rest of the cycle, the 
latge vortices are progressively washed away and the flow reattaches on the 
upper surface. The center of pressure still has large excursions: this would 
be felt as buffeting in an airplane. 

Figure 25 presents the numerical results for the normal force coefficient and 
the moment coefficient, compared with experimental results reported in [60]. 
The experimental results are phase werrged over a large number of cycles. 
Three successive cycles of computation are shown; the first cycle is not as 
representative since it smarted from fully irrotational flow. Three cycles are 
too small a sample for phase-averaging and this is why individual cycles are 
shown. 

The agreement is quite good. However, the peak value of the normal 
force is lower in the computations. The peak value of the moment is also 
significantly lower; the experic ents exhibit a large peak during the phase 
called "moment stair. Another area of disagreement is during the low 
incidence period. Reattachment often seems to take an unexpectedly long 
time 5 n the computations and this might account for the difference. It is 
hoped tnat these disagreements will be better understood when the tunnel 
walls 2 re included in the computation. 

Figure 26 slices the pressure distribution on the airfoil during dynamic 
stall, after it has been Fourier-transformed in time. The mean pressure is 
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shown, as well as the in-phase component (the phase being given for the 
incidence by Eq. 99) and the out-of-phase component. Each component has 
been Don-dimensionalized, the mean is divided by the mean incidence ao 
and the first harmonics are divided by the amplitude a*. The agreement 
between data extracted from different cycles is fair. The difference between 
the mean and the in-phase components, and the existence of an out-of- phase 
component, are manifestations of the non-linear behavior of the flow. 

d) Tilt-Rotor wing in hover. 

The "Tilt-Rotor” concept combines many of the advantages of wingborne 
aircraft and of helicopters and has been studied for many years. NASA and 
Beil have recently conducted a very successful experimental program which 
included the construction and testing of the XVI 5 aircraft. A production 
version called JVX is now being developed and more emphasis is being placed 
on performance. One factor in the hover performance is the download ex- 
perienced by the wing in the downwash of the rotors; both an experimental 
program and a numerical study are under way at NASA Ames to improve 
this aspect of the system. 

A Tilt-Rotor aircraft has two rotors that can* point vertically for vertical, 
helicopter-like take-ofT and landing, and also point forward for airplane-like 
forward flight. The Tilt-Rotor is much more efficient than a helicopter in for- 
ward flight because it avoids the asymmetric conditions than deteriorate the 
operation of the helicopter rotor, with the advancing blade experiencing drag 
rise from high subsonic Mach numbers and the retreating blade experiencing 
low dynamic pressures and stall. As a result the Tilt-Rotor is much faster 
and has a lower fuel consumption than the helicopter; its development was 
delayed mostly by aeroelasticity and control problems. 

The Tilt- Wing concept, in which the wing tilts with the rotors, was not 
as successful probably because it required a heavy articulation at the wing 
root, and had more severe control problems. However, for a Tilt-Rotor when 
the .otors are pointing up the wing is still horizontal and placed in the 
downwash of the rotors; this creates a downward force on the wing which is 
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very significant especially since hover is the most critical condition in terms 
of payload. The objective is to decrease this force by improving the shape of 
the clean wing and by adding devices, like flaps, that could also serve as high 
lift devices for slow forward flight. 

Almost all of the wing is in the downwash, which is far from being uniform. 
An accurate computation of the three dimensional flow is presently out of the 
question; the computations will be done in two dimensions and are expected 
to show trends and provide guidelines for the choice of the airfoil and landing 
devices on the wing. 

The Tilt-Rotor wing can be quite thick since the Mach numbers it reaches 
are not very high; the main advantage is a lighter structure and more space for 
fuel and other components, but the larger thickness coincidentally improves 
the situation with regard to download. The NACA 4421 airfoil was chosen 
for the first stage of the numerical study and three landing devices were 
considered: two kinds of trailing edge flap and a leading edge device. Both 
flaps span 25% of the chord; the first one pivots about a hinge that is within 
the thickness of the airfoil and the second one has the hinge at the lower 
surface. The deflection of the flap introduces an arc of a circle as part of the 
top surface; the first flap has a circle of smaller radius than the second. This 
curvature was expected to strongly influence the separation of the boundary 
layer and therefore the global force. The leading edge modification involved 
removal of about 10% of the chord, thus reducing the effective chord and 
increasing the leading edge radius, with again a favorable effect on separation 
(the part of the airfoil that is removed would be placed under the leading edge 
and thus shielded from the mainstream). 

The KPD2 program was chosen. All cases were run with the same 
parameters to allow a comparison. The Reynolds number was 10 7 , and the 
boundary layer underwent transition; all the shapes were defined by about 
200 points, with At = 0.03, R 0 = c = 0.02, Dq = 0.2 and N v = 1000. 
Several thousand time steps were necessary for an average value to emerge; 
this represents 20 to 30 minutes on the CRAY-1. 

The flat plate was treated under the same conditions to provide a reference; 
the result for the drag was about 3.5, compared to the experimental value of 
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about 2 [5]. Thus all the results are probably higher than the correct values; 
however, this does not prevent a comparison between different configurations, 
which was the main objective of the study. 

Figure 27 shows the drag computed for the flat piate reference, both flaps 
with deflections ranging from 0 to 9(f, and the modified leading edge with 
the large radius flap. Tb f clean airfoil had a drag coefficient of 2.8, which is 
also thought to be higher than the correct value. 

It appears that the two flaps give results that are much closer than one 
would have expected, and that the small radium flap can even have less drag. 
In retrospect, the reason is probably that even when the flow separates from 
the top surface at the hinge it can reattach on the flap; see Fig. 28 (in 
Fig. 28 the airfoil, vortices and force are shown, as well as the instantaneous 
streamlines). The size of the separation bubble will not influence the drag. 
Furthermore, at low deflections the large radius flap actually increases the 
chord of the airfoil, compared to the small radius flap; this could explain why 
the small radius flap is better below 45°. At higher deflections the advantage 
of the large radius flap in terms of delaying the separation seems to come 
into play; however the difference is probably not large enough to override the 
structural considerations and dictate the choice of the flap. 

Another lesson learned is that low flap deflections barely reduce the drag; 
here, a linear decrea»e could have been expected. The drag then reaches 
a minimum around 75°, and then rises again; this is probably because the 
flow does not reattach on the flap any more (Fig. 28). This character of 
the drag curve, with a plateau at low deflections and a minimum at less 
than 90 P, was observed in measurements performed on the aircraft itself; a 
more quantitative comparison will not be attempted since the flow is strongly 
three-dimensional . 

A third lesson learned is that the leading edge modification reduces the 
drag appreciably without flap deflection, but loses almost all of its effect 
when the flap is deflected to 75° (Fig. 29). This was disappointing and no 
convincing explanation has been found. Other types of leading edge devices 
are now considered. 

The flows exhibited a shedding of large vortices; the Strouhal numbers were 
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of the order of 0.15, which represents rather slow variations of the loads. Some 
configurations also gave rise to a significant force normal to the airflow. The 
clean airfoil had low pressure near the leading edge and a small forward force. 
With the flap deflected, rearward forces of the order of 0.5 were observed. 
These might have to be taken into account since the available control power 
is not very high in a hover condition. 

As a whole this study demonstrated the flexibility of the method in treating 
many different shapes, and produced consistent and useful results even in a 
"first pass” with a set of parameters that was not optimal. It also gave clear 
indications about the structure of the flow, for instance the presence of a 
bubble at the hinge. In the near future a more extensive study will involve 
the actual XV15 airfoil, flaps and other possible devices, wall effects, and will 
be directly compared to wind tunnel tests. 


5) Results from the program KPD3. 

The program KPDZ has been applied only to the circular cylinder. 
Attempts to treat airfoils with KPDS were unsuccessful; the algorithm seems 
to be unable to treat the trailing edge region properly. The most likely ex- 
planation is that the boundary layer assumptions are simply not applicable 
in a region where the wall has so much curvature. 

The flow around a circular cylinder is the most classical bluff body problem 
and the best documented. It is also a good test of the ability of the program 
to predict the stucture of the boundary layer. Its dependence on the Reynolds 
number is known to be strong. The range of Reynolds numbers that were 
considered is from 10 4 to 3.16 x 10 7 . 

The value 10 4 is considered to be roughly the lower limit because at Re — 
10 4 the thickness of the viscous region is about 0.035, which is not very 
small any more compared to the radius of cutvaturo. which is 1. As the 
Reynolds number decreases, use of the boundary lay?r assumptions become 
less justifiable. 

The value of 3.16 x 10 7 was chosen as the upper limit; experiments nave 
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been conducted up to about 0.8 x IQ 7 . Agreement between numerical and 
experimental results was found to be quite good in the upper range. Therefore 
it seems reasonable to consider extending the range of the numerical simula- 
tions by this factor of 4 in Reynolds number. Furthermore, the turbulence 
model has been validated at Reynolds numbers slightly higher than 10 7 [50]. 

Fig. 30 shows the global numerical and experimental results for the drag of 
a circular cylinder. Fig. 31 shows the Strouhal number and Fig. 32 the mean 
separation angle. Fig. 33 shows the pressure distribution at eight different 
Reynolds numbers. These cases include Re = 10 4 , 10 5 , 10 8 and 10 7 , some 
extra cases in the critical region: 10 5,5 (3.16 x 10 5 ), 10 5,75 (5.62 X 10 s ) and 
10 6 - 5 (3.16 x 10 8 ), and Anally 10 7,5 (3.16 x 10 7 ). Experimental results are also 
shown when available, [61], [62], [63], [64]. Fig. 34 shows the wall shear stress 
coefficient at the same Reynolds numbers. In order to make comparisons 
easier the scales are the same for all the pressure plots and for all the shear 
stress plots, except 10 4 . 

The Strouhal number was determined by counting the apparent number of 
periods in the lift signal during the simulation, using plots like the ones in 
Fig. 40. II the length of the sample is of the order of 5 periods, this way of 
determining the Strouhal number is accurate although somewhat arbitrary. 
Another way is to Fourier-analyze the signal; however one still has to " clip” 
the sample to simulate a periodic signal, which is also arbitrary. The Fourier 
analysis method is reliable only if long samples are available; experimental 
results often provide hundreds of cycles, but computations do not, for obvious 
cost reason?. As a result, the spectra contain subharmonics which would 
probably disappear if longer samples were available. In that sense, these 
harmonics are not meaningful, physically (see Fig. 41). 

The mean separation angle was defined as the time-average of the angle at 
which the instantaneous shear stress is zero. This does not coincide with the 
angle at which the time-average of the shear stress is zero, because the shear 
stress is a strongly non-linear function in that region. Actually, in some cases 
the numerical results for the shear stress do not cross zero. This point will 
be developed later. 
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a) Subcritical regime. 

This regime exists for Reynolds numbers under approximatively 10 s . The 
drag and Strouhal number are almost constant, and the evolution of the flow 
as the Reynolds number increases is weak (although the base pressure does 
change, and the lack of change in the drag is fortuitous [65]). Both boundary 
layers separate while still laminar and the separation point does not move 
significantly. The shear layers transition in the wake and transition occurs 
closer to the cylinder as the Reynolds number increases [65]. Apparently the 
"drag crisis” starts around Re = 2 x 10 s when the process is complete and 
transition occurs in the boundary layer just before or after separation, rather 
than in the wake. As a result, the boundary layers reattach and separate 
again, only farther downstream. 

In the numerical results the boundary layers never transitioned for 
Reynolds numbers of less than 3.16 x 1C\ Thus the only evolution with 
Reynolds number is the thickness 6 of the viscous region, which scales with 
Re”" 1 / 2 , from 0.035 at 10 4 to 0.015 around 10 s . The average separation 
angle is almost constant around 84°. The Strouhal number is between 0.205 
and 0.195, which is considered a good prediction. The drag coefficient on 
the other hand slowly decreases from 1.42 at 10 4 to 1.03 at 10 s and 0.88 at 
3.16 x 10 5 . This is not very good since, experimentally, the drag does not 
decrease in that range. 

The pressure and shear stress distributions are correct qualitatively, but 
the quantitative differences which cause the inaccuracy in the predicted dr~g 
are obvious. At 10 s the bare pressm is -0.85, compared to -1.25 in the 
experiments. There is a significant pressure rise downstream of the separation 
point, which is not correct physically. 

The average shear stress remains positive until the 11CP point is reached, 
while the mean separation angle is 85°. The reason for this paradox is 
clear on Fig. 35b, which shows a still of the simulation at Re = 10 s . 
The instantaneous shear stress crosses zero, but has significant oscillations 
downstream of that point. As the flow oscillates the separation point moves 
back and forth and, in the average, the shear stress remains very close to 
zero instead of being frankly negative as in the experiments. This value 
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of the shear stress is small and results from the averaging of much larger 
numbers; therefore, although the situation seems inconsistent physically, the 
inaccuracy is quantitatively small. 

Figure 35a shows the state of the inner flow at the same instant. It includes 
vorticity contours and the values of U e , V e and C p (the inner region has been 
unwrapped from the cylinder and the scale expanded in the n direction). 

Starting from the center of the figure (front of the cylinder) the two boun- 
dary layers of opposite sign divide and progress along the walls. The mag- 
nitude of U e increases, the pressure decreases (favorable gradient) and V e is 
negative (flow into the inner region). The wall shear stress (plotted in Fig. 
35b) reaches a maximum . In the region near 90° the situation changes. The 
pressure gradient reverses and is now adverse; V t becomes positive and the 
flow as a whole starts moving away from the wall, taking the vorticity with 
it. The shear stress rapidly falls to zero. Finally the vorticity leaves the 
inner region: XJ e falls to zero and the shear layer is now in the outer region 
under the form of vortices. Downstream of the separation region the pressure 
fluctuates as the eddies contained in the outer region progress along the wall. 
However the time-averaged pressure is smooth and almost constant in the 
wake region. 

The behavior that was just described is consistent with what is known of 
the flow. The comparison between Fig. 35b and Fig. 36 probably explains 
why the numerical results differ significantly between Re = 10 4 and 10 5 . 
Fig. 36 shows that at 10 1 the separating shear layer is much thicker than at 
10 5 . As a result, its breakdown into circular vortices is slower and occurs 
farther away from the wall region. Being thicker the shear layer is also 
better represented by the limited number of vortices that are available. In 
contrast, at 10 5 it becomes difficult for the vortices to properly model the 
shear layer; the vortices are seen to linger near the wall while they separate 
more cleanly at 10 4 . Very probably, the vortices cause too much mixing 
too soon, which creates a situation analogous to a turbulent boundary layer: 
a layer with strong mixing and a laminar sublayer. This situation is not 
incorrect physically, but the intense mixing is partly of numerical origin and 
may not be quantitatively correct. 


81 



To support this conjecture the flow at Rt' = 10 5 was simulated with two 
levels of numerical resolution. At the low level N v = 1000, N w — 200 and 
Ai = 0.025 were used. At the high level N v = 1600, N w — 300 and A t = 
0.02 were used. The results were significantly different: the Cd is 0.8 with low 
resolution and 1.03 with high resolution (the correct value is 1.2). We can 
conclude that the regime just below critical and the critical regime itself are 
the most difficult to simulate, because the shear layers are very thin, and the 
level of resolution used is not quite sufficient. 

Other quantities have been computed to compare with experiments. The 
average velocity at the edge of the attached boundary layer is shown in Fig. 
37. The agreement with the experiments reported in [65] is excellent. This 
is consistent with the good agreement shown by the pressure in the same 
region. 

The value of the streamwise velocity on the centerline behind the cylinder is 
shown on Fig. 38 and compared with experiments from [62]. The agreement 
is good in the near wake but worsens in the far wake. The computations 
seem to introduce less dissipation than the experiments indicate. This was to 
be expected since tne successive merging of the vortices in the wake tends to 
concentrate the vorticity whereas in the real flow the turbulent stresses spread 
it. It seems that the description of the wake is adequate up to approximatively 
2 diameters of the back of the cylinder. 

In conclusion, the subcritical regime, even though the boundary layers are 
laminar, appears to be more sensitive to transition in the near wake thaD 
was expected. A more accurate description of this region is probably needed. 
Another possible source of error is the three-dimensional character of the 
real flow; the consensus seems to be that large scale three-dimensional effects 
are not very strong in the subcritical regime, but small scale effects like 
streamwise vortices might play an important role in the transition region. 

b) Critical regime. 

This regime is the most complex and the most difficult to model. An 
intense experimental activity to measure and describe the flow accurately is 
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being carried out [61]. It has long been known that tue drag coefficient takes 
very low values, of the order of 0.25, and that vortex shedding is disrupted 
[64]. It was also known that this "drag crisis” is associated with transition in 
the boundary layers [1]. Recent experiments have shown in addition that the 
decrease in the drag is probably not smooth but occurs in steps instead, and 
that the flow can be in states with a non-zero average lift [61]. Finally, the 
span-wise coherence of the real, three-dimensional flow is probably weakest 
in the critical regime [67]. 

The numerical method essentially failed to model the critical regime. This 
is not very surprising; the dramatic changes in the flow pattern are probably 
associated to reattachment of the boundary layers, or of only one of them. 
Reattachment is a delicate phenomenon to simulate; it requires a very ac- 
curate coupling if a zonal viscous-inviscid approach is used, and in any case 
a very fine transition and turbulence model [55]. The turbulence model used 
here was not designed with low Reynolds numbers in mind. Furthermore 
the transition model was designed for attached boundary layers. It is a "one- 
dimensional” model in which the whole boundary layer transitions at once. 
In contrast the results in [55] showed a very non-uniform turbulent energy 
across the layer. 

As mentioned earlier transition was not predicted by the simuhMons for 
Reynolds numbers of 3.16 x iO 5 or less; for Reynolds numbers of 10 8 or more 
transition always occurred in both boundary layers. In that sense, the critical 
regime predicted by the computations is between these two values of Re. At 
Re s= 5.62 X 10 s transition was intermittent; in Fig. 39 transition is seen to 
occur in only one of the boundary layers, as can be inferred by the sudden 
rise in the wall shear stress. However the simulation did not stay locked in 
either position (Fig. 39a and 39b). As a result the stable, asymmetric, lifting 
situation that was observed in experiments [61] was not obtained. The low 
drag values were not observed either. 

'i he shear stress distribution at 5.62 X 10 5 reflects this intermittent tran- 
sition. it has a plateau between 90° and HOP . The pressure distribution does 
not compare well with experiments, if only because it is essentially symmetri- 
cal. 
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One feature of the critical regime that is observed is the alteration of 
the vortex shedding pattern. Fig. 40 shows the drag and Uft as function 
of time in the simulations at Re = 10 4 , 3.16 x 10 5 and 3.16 x 10 6 , and 
Fig.41 shows the spectra of the lift signals. The lift has a fail iy well defined 
oscillation at 10 4 (although strong modulations are evident) and very regular 
oscillations at 3.16 x 10 6 . At 3.16 X 10 5 the lift keens the same sign for 
much longer periods of time and its behavior is very far from being harmonic. 
This impression is confirmed by the examination of the spectra: the peak 
is much broader at Re — 3.16 X 10 s . This might be an indication that a 
small modification of the conditions could cause the simulation to find the 
stable asymmetric configuration. However, the samples were too short for the 
spectra to be entirely reliable. For instance, the peak at a Strouhal number 
of approximatively 0.08, at Re = 10 5 , is thought to be spurious. 

In conclusion, the correct simulation of the critical regime would probably 
require a more elaborate transition and turbulence model (the McDonald-Fish 
model used for the simulations in [55] was quite complex). However it seems 
that the method should first be improved until it simulates the subcritical 
rejme very accurately before another attempt is made on the critical regime. 

c) Supercritical regime. 

This regime exists for Reynolds numbers above approximatively 4 x 10 6 . 
Both boundary layers transition before separating and the flow is similar to 
the subcritical flow, except that separation occurs much later, whrh results 
in a narrower wake, a lower drag coefficient and a higher shedding frequency. 
The drag coefficient seems to be almost constant again in that range (although 
the experiments do not exceed 8 x 10 6 , so that the "plateau" is quite narrow, 
covering oniy a factor 2 in Reynolds number). 

The numerical results show slow variations of drag coefficient, Strouhal 
number and separation angle up to Re = 3.16 x 10 7 . Separation does 
occur later than in the subcritical regime, and the changes in wake pattern, 
shedding frequency and drag are all correct qualitatively. The pressure and 
friction distributions are almost Reynolds number-independent; the transition 



point just moves upstream as the Reynolds number increases. The agreement 
•with experiments for the pressure distribution is not excellent and the drag 
coefficient is lower than the experiments indicate. However there is significant 
scatter in the experimental results. The agreement in the shetr stress is poor; 
the experiments do not show the steep rise associated with transition. 

Numerically the simulations look "healthier” in that range than at lower 
Reynolds numbers. The separation is frank and the pressure distribution is 
smooth (Fig. 42 shows a still of the inner flow and Fig. 43 four stills of 
the outer flow. The non-dimensional time, based on velocity and radius, is 
indicated and the four stills cover approximative^ one period of shedding). 
The reason is probably that the separating shear layer is thicker again, due to 
the Reynolds stresses in the inner region. Another sign is that at Re = 10 6 
the simulations with low and high resolution agree very well (Figs. 33 and 
34). This indicates that the resolution is sufficient. 

In conclusion, the supercritical regime is easier to simulate from a numerical 
point of view, and the Baldwin-Lomax turbulence model seems to work weli. 
The accuracy of the drag is not easy to assess since few experiments have 
been conducted in that range. 

This concludes the description of the result? obtained from the programs 
KPD1, KPD2 and KPD3. These results show that the method can 
reproduce many of the features of separated flows a' ' is generally in fair 
to good agreement with experimental predictions. There is, however, room 
for improvement of the quantitative agreement, even in the framework of a 
two-dimensional method. 
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V) CONCLUSION. 


In this work we have developed a new numerical method for two- 
dimensional separated flows. It incorporates an improved version of the 
Vortex Method, to treat the inviscid outer region, and widely used integral 
or finite difference methods, to treat the viscous inner region. 

The numerical method uses a formulation of the equations in terms of the 
vorticity. This formulation has been shown to be mathematically equivalent 
to the conventional velocity/pressure formulation and includes cases with 
solids in non-uniform motion or rotation. It could also treat bodies in a 
uniform shear flow. 

The flow is treated as inviscid away from the solid walls. This is motived 
both by physical arguments (in the wake the large structures dominate and 
are not very sensitive to the viscosity) and by numerical arguments (it is 
not practical to model structures so fine that the viscosity influences them 
significantly). 

It is shown that the Vortex Method can produce a significant numerical 
diffusion. This diffusion is present only if velocity gradients are present; this 
is an advantage, compared to E’llerian methods. However strong gradients 
are present in most viscous flows, especially near the wall. As a result, the 
numerical diffusion can have a dominant influence on the simulation. 

In the first main version of the algorithm all the vorticity is represented 
by vortices, but their departure from the wall is delayed, 1/ necessary, so 
that separation occurs at the proper location. This location is predicted by a 
boundary layer solution based on an integral method, which is run in parallel 
with the vortex solution. 
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In the second main version, a thin region near the wall is treated ^.s viscous, 
with boundary layer assumptions. This inner region is treated by a Finite 
Difference method. The key assumption is not the neglect of the streamwise 
viscous term ud 2 ui/ds 2 , since centered finite differences are used in the s 
direction anyhow, it is the neglect of the dv/dn term in t v e vortic»ty (to avoid 
the solution of an elliptic equatiou for the velocity) and the coupling with tke 
outer flow by means of a vortex sheet (to simplify the interaction between 
the two regions). The coupling algorithm is the most delicate element of the 
algorithm, and ca^ neate problems especially if the wall has tight curves. 
However it is an essential ingredient to make the method versatile and lessen 
its dependence on empiricism. 

The new method conserves the traditional advantages of the Vortex Method 
for the treatment of the wake: it is grid free and accurate in modeling 
transport phenomena, it treats the far field in a simple and accurate wav 
It requii j much le > memory space and possibly less computing effort than 
comparable Finite Difference methods . In addition the method can now 
treat arbitrary solid bodies or grou is of bodies; conformal mappings are not 
involved. 

Many choices have to be made when designing a practical method. Some 
of these choices could be made in a systematic way: for instance the design 
of the merging device. Some choices were made on a more intuitive basis, 
and therefore could probably be improved: for instance some details of the 
coupling procedure, and of the implementation of the turbulence model. 

In its first version the new method treats bluff bodies reliably and is quite 
accurate. With separation properly controlled it simulates th* stall of air- 
foils accurately; the comparison with experimental results is very encourag- 
ing. Some possible reasons for ,he remaining disagreement are the thrt a - 
dimensional character of the real flow, which the method is unable to account 
for, and the interference with wind-tunnei walls. These wall effects are being 
added to the method. 
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The second version has been used for an extensive study of the flow past 
a circular cylinder at Reynolds numbers ranging from 10 4 to 3. x 10 7 . The 
change in wake pattern and decrease in the drag coefficient as the boundary 
layers switch from laminar to turbulent are clearly observed, even though 
the pattern in the critical range itself is not correct. In this range, the 
+ Iiree-dimer^ional effects are probably quite strong, and the modeling of 
transition and turbulence is the most delicate. In the supercritical regime 
( - oove 3. x 10 6 ) most of the flow characteristics are essentially constant and 
agree well with experimental results. In the subcritical regime (from 10 4 to 
2. x 10 5 ) the shedding frequency is accurate and the drag is close to the 
correct value. However the drag coefficient tends to decrease steadily as the 
Reynolds number increases, which is not correct. This seems to be due to a 
"transition" of the separating shear layers which is not accurate, physically. 
This problem could be alleviated significantly, although not suppressed, by 
refining the discretisation. However, a further increase in the computation 
cost might not be the best answer; it would be preferable to improve the 
algorithm in some way. 

The prospects for the method developed here to be applied to practical 
cases, at least in its first version, are good. This version already has two active 
research applications: the Tilt-Rotor Wing study and the Vortex Flowmeter 
study. Both projects are continuing and a favorable interaction between the 
numerical work and the experimental work is building up. The second version 
will probably be used for further investigation of the circular cylinder flow. 

For the future many extensions are possible, especially for the second 
method. They include: 

• The introduction of a boundary layer model of higher order, with curva- 
ture terms, in the hope of treating bodies with tight curves, 

• The refinement of the coupling algorithm in the spirit of the " strong in- 
teraction” theories for the separation region. A more mathematical approach 
could be taken; the improvement of the coupling might involve a more com- 
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plex model than the single vortex sheet placed at 8 *. 

• In the long term the inner region could be of order 1 in thickness and be 
treated by a full Navier Stokes solver. This should remove all the problems 
associated with boundary layer assumptions. On the other hand it would 
probably demand a much more complex coupling procedure, as well as the 
solution of an elliptic equation in the i^ner region instead of Eqs. 29 and 30. 

• The introduction of a more elaborate model for transition and for the 
turbulent stresses, in the hope of reproducing the well-known "drag crisis" of 
the circular cylinder and the recent findings about stable asymmetric states 
and discontinuous variation of the drag coefficient. However, such a model 
is not readily available, and in addition it might be that a two-dimensional 
model will never be able to account for the drag crisis. This question is open. 

• It has been proposed to couple the Vortex Method to a finite difference 
method to treat compressible flows; the finite differences would handle the 
compressible effects efficiently on a fairly coarse mesh while the vortices would 
provide a fine description of the thin shear layers in the wake. The theoretical 
work to support these ideas has not beei done yet. 

• The extension to three dimensions, to treat wings or cars for instance. 
The validity of two-dimensional analysis for such cases is very lirui'ed, and a 
relatively crude three-dimensional method might be more useful than a ” fine- 
tuned” two dimensional method. The power of a CRAY should be sufficient 
for simple cases, and the treatment of the solid should not be very different 
from the two-dimensional case. What is needed is a new three-dimensional 
Vortex discretization that would be reliable and more flexible than the present 
methods, with their connected filaments; such methods are currently being 
developed. 


APPENDIX A 


This appendix contains the description of the versions KPDl and KPD2 
of the program. They differ from KPD3 in the way they handle the wall 
region. 

1) Program KPDl. 

This version is a pure Vortex Method. It is mentioned mostly to clarify 
some issues about the Vortex Method and its ability to simulate some viscous 
flows while solving the inviscid equation. 

In KPDl even the boundary layer vorticity is represented by vortices that 
belong to the outer flow. The vorticity is immediately injected into the vortex 
region instead of transiting through an inner region as in KPD3; the new 
vortices follow the wall until they separate. The treatment of the wall region 
is thus very simple (Fig. 36). The rest of the algorithm (blob shape, time 
integration, merging, etc.) is the same as in the ou er region of KPD3 (see 
Chapter IV). 

Fig. 37 is a flow chart of KPDl; it is quite a simple program and a listi^ 
of it is included in Appendix B. 

2) Program KPDl. 

In this version, the premature separation of the vortices that often occurs 
with KPDl (see Results), is artificially prevented on the basis of a boundary 
layer computation. The boundary layer is compute by an integral method. 

The integral method is chosen for its simplicity and low cost. The solvers 
used are taken from Ref. [68]: Thwaites’ nr jiod for the laminar part of the 
boundary layer and Head’s method for r .e turbulent part. Both methods 
use a small number of integral quantities in the boundary layer as degrees of 
freedom and compute the boundary layer by marching in the stream direction, 
using empirical laws for tb* .volution of the integral quantities. Neither 
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method can be used downstream of the separation point. Both methods were 
designed for steady boundary layers; the unsteady effects on the boundary 
layer were neglected. Finer integral methods exist that include the unsteady 
effects; however these are probably weak compared to the interaction with 
the outer flow vhen massive separation occurs. Transition was modeled by 
the Sclicht' g-Granville criterion as in KPD3. 

The sapling between the vortex flow and the boundary layer is done as 
folic;? s: the vortex flow determines the instantaneous pressure distribution, 
?nd the boundary layer determines the separation points. Fig. 38 is a flow 
chart of KPD2. 

At each time step, the pressure distribution is computed from the outer 
flow solution, using Equation 13. This pressure is then used as the forcing 
function to solve the boundary layer equation (this is upstream of separation 
and Bernoulli’s theorem applies). The upper and lower boundary layers are 
computed, from the attachment point to their first separation point. 

Experience with KPDl shows that the vortices, if they are left free, always 
separate too soon. Therefore making them separate it the proper location 
is simply a matter of delaying their natural separation. To achieve this, all 
the vortices that were created upstream of the desired separation point are 
marked as ’’temporary” and after one time step are removed and replaced 
by a fresh layer of vortices. This prevents the vorticity layer from thickening 
prematurely. The vortices that are created downstream of the separation 
point are marked as "permanent” and treated as in KPDl ; being in an 
adverse pressure gradient, near the separation point, they leave the vicinity 
of the wall quite promptly. Thus the vorticity layer is effectively released 
into the outer flow a short distance downstream of the separation that was 
predicted by the boundary layer solver. 
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PROGRAM KPD1(1NPUT,0UTPU i , TAPE 1.TAPE8, TAPED, 

> TAPE5=INPUT, TAPE6=OUTPUT) 

C 

C 2D VORTEX TRACING SIMULATION PROGRAM. 

C 

C PROGRAM WRITTEN BY PHILIPPE R. SPALART 
C 

C IN COLLABORATION WITH ANTHONY LEONARD. 

C 

C SEE AIAA PAPER 81-1246. 

C 

C NASA AMES RESEARCH CENTER, NOVEMBER 1982. 

C 

C INCOMPRESSIBLE FLUID. UNIFORM VELOCITY U1NF 
C AT INFINITY, UINF IS COMPLEX (MAGNITUDE: ABSU1N, 

C INCIDENCE IN DEGREES: ALPHA). 

C THE PROGRAM COMPUTES THE UNSTEADY FLOW, 

C STARTING FROM POTENTIAL FLOW. 

C THE SOLID SHAPE IS ARBITRARY, GIVEN BY THE ROUTINE SOLID. 
C IT CAN BE MADE OF SEVERAL SEPARATE BODIES. 

C THIS VERSION OF THE PROGRAM IS SIMPLIFIED 
C AND MORE ROBUST (COMPARED TO THE ONE USED 
C FORTHE PAPER) AND DOES NOT HAVE ANT' 

C REYNOLDS NUMBER EFFECTS. 

C THEREFORE IT IS SUITED FOR SHAPES THAT DO NOT 
C HAVE A STRONG REYNOLDS NUMBER DEPENDENCE, 

C ESPECIALLY SHAPES WITH SHARP CORNERS. 

C FOR OTHER SHAPES, CIRCLES FOR EXAMPLE, 

C THIS PROGRAM GENERALLY GIVES RESULTS 
C CLOSE TO LAMINAR RESULTS (RE=10**5 OR SO). 

C 

C THIS MAIN PROGRAM CALLS 
C READPR, GMTRY, INIT, MERGE, BCBODY .MOVE. 

C WRITING RESULTS ON TAPE8, THEN READING THEM 
C FROM TAFE9 ALLOWS RUNS TO CONTINUE EACH OTHER. 

C 

COMPLEX Z,V,FORC£(2), HUB, ZO, WALL, ZCR.VM, UINF AVFO 
REAL MOM(2),DPDS(215,l),CP(2l5,l),CPAV(215) 

C 

COMMON/VORTEX/NVORT,Z(2000),V(2000),VM(20()0), 

> GAMMA(2000) 

COMMON/SOLID/NBDIES,NWALL(2),WALL(2l5,l), 

> ZCR(215,1),HUB(2),Z0(2),THETA(215,1),NINC(2) 

> JNC(15,2),XMAX(2) 
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ORIGINAL FALSE ** 
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COMMON/SYSTEM/ND IM,A(215, 215) ,X(2 15), IP VT(215) 

READ AND PRINT THE PARAMETERS. 

CALL READPR(ISTART,NDES,NSTEP,N2,ABSUIN, ALPHA, 

> DELT ,D0,GAMMA0) 

SET UP THE GEOMETRY 

DEFINE THE SOLID AND THE CREATION POINTS ON ITS SURFACE. 
COMPUTE AND GAUSS ELIMINATE MATRIX OF INFLUENCE 
COEFFICIENTS BETWEEN WALL POINTS, AND DO OTHER 
THINGS THAT DEPEND ONLY ON THE SOLID. 

CALL GMTRY(SIGMA2, CHARD) 

INITIALIZE THE TIME-DEPENDENT VARIABLES. 

CALL INIT(ISTART,NSTART,T,VO,ABSUIN, CHARD, DO) 

UINF=ABSUIN*CEXP(CMPLX(0 ., ALPHA* ATAN( 1 .) /45 .)) 

AVFO=0 


MAIN LOOP; ADVANCE FLOW TIME STEP BY TIME STEP. 
WRITE (6,104) 

I FORMAT(//, B STEP BY STEP EVOLUTION OF THE FLOW:",//) 

NEND=NSTART+NSTEP-1 
DOl N=NSTART ; NEND 

MERGE VORTICES TO KEEP THEIR NUMBER REASONABLE. 
CALL MERGE(D0,V0,NDES,N) 

TREAT BOUNDARY CONDITION AT THE BODY BY AN 
EXCHANGE OF VORTICITY. 

CALL BCBODY(FORCE,MOM,GAMMAO,N,T,UINF, 

> CHARD, NOLD, CP, DELT .SIGMA2) 

AVFO=AVFO +FORCE( 1 ) 

MOVE VORTICES. 

CALL MOVE(UI NF, SI GMA2, NOLD, DELT) 

END OF MAIN LOOP. 

STORE RESULTS IN CASE WE WANT A FOLLOW UP TO 
THIS RUN. 

WRITE(8) N,T,NVORT,Z, GAMMA, VM, VO 
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C OUTPUT AVERAGE LOADS. 

AVFO=AVFO/NSTEP 

WRITE (G,100)REAL(AVFO)^IMAG(AVFO) 

100 FORMAT (/,24H AVERAGE DRAG AND LIFT: ,2F8.4) 

C 

C END OF RUN. 

STOP 

END 

SUBROUTINE READPR(ISTART,NDES,NSTEP,N2, 

> ABSUINALPHA,DELT,DO,GAMMAO) 

C 

C READ PARAMETERS. 

C ISTART=I IF RUN IS FROM TIME 0. 

C IS TART =0 IF IT IS A FOLLOW-UP. 

C NDES DESIRED NUMBER OF VORTICES. 

C THE PROGRAM WILL ROUGHLY MAINTAIN 

C THE NUMBER OF VORTICES AT NDES. NSTEP NUMBER OF STEPS. 

C ABSUIN MODULUS OF UINF, ALPHA INCIDENCE 
C IN DEGREES. DELT TIME STEP. 

C DO PARAMETER IN MERGING DEVICE. DO SMALLER PUTS MORE 
C VORTICES NEAR THE SOLID AND LESS FAR FROM IT. 

C GAMMAO ALLOWS THE USER TO DISTURB THE FLOW TO 
C MAKE IT REACH THE SHEDDING REGIME FASTER. 

C GAMMA0=0 LEAVES IT UNDISTURBED. 

C GAMMAO. NE.O ARTIFICIALLY ADDS A CIRCULATION 
C GAMMAO AT THE BEGINNING OF THE RUN. 

C (GAMMAO IS IGNORED IF ISTART=0). 

READ (5,200)ISTART,NDES,NSTEP,N2 

200 FORMAT (1 1,315) 

READ (5, 201)ABSUIN ALPHA, DELT ,D0,GAMMA0 

201 FORMAT (5F8.5) 

C 

C PRINT THE PARAMETERS. 

WRITE(6,100) 

100 FORMAT(//," VORTEX SIMULATION OF BLUFF BODY FLOW” ,//) 
IF(ISTART.NE.0)WR1TE(6,101)GAMMA0 

101 FORMAT (" THIS RUN STARTED WITH CIRCULATION: ",E8.2) 
WRITE (6,102)NDES 

102 FORMAT(/," APPROXIMATE NUMBER OF VORTICES: ” ,16) 

WRITE (6,103)ABSUINALPHA 

103 FORMAT (/, 

>" FRLESTREAM VELOCITY MAGNITUDE AND INCIDENCE ”,2F7.4) 
WRITE (6,111)DELT 
111 FORMAT(/,llH TIME STEP ,F7.4) 


05 



o o no 


ORIGINAL PACE IS 
OF POOR QUALITY 


WRITE (6,H3)D0 
113 FORMAT (/, 

>" CHARACTERISTIC DIMENSION IN MERGING DEVICE ",F7.4) 

C 

RETURN 

END 

SUBROUTINE GMTRY(SIGMA2,CHARD,LA,I<A,LB,KB) 

C 

C OBTAIN SOLID SHAPE, COMPUTE SOLID-RELATED ARRAYS, 

C GAUSS ELIMINATE THE MATRIX, -p MISCELLANEOUS. 

C NBDIES NUMBER OF BODIES. 

C NWALL NUMBER OF WALL POINTS ON EACH OF THEM. 

C WALL ARRAY OF WALL POINTS. ZCR ARRAY OF CREATION POINTS. 
C THETA POLAR ANGLE OF ZCR POINTS. 

C THETA IS USED TO FIND IF VERTICES ARE 
C INSIDE SOLID. ZO USED ALONG WITH THETA. 

C HUB(L) IS THE HUB OF THE BODY "L” . INC WILL 
C HELP FIND OVER WHICH WALL POINT THE VORTEX IS, 

C BY BISECTION. THE FIRST DIMENSION 
C OF INC MUST BE AT LEAST LOG2(NWALL). 

C A IS THE MATRIX OF INFLUENCE COEFFICIENTS FROM 
C CREATION POINTS TO WALL POINTS. 

C 

INTEGER FIRST, NEXT(215) 

COMPLEX WALL,ZCR,ZO,HUB,INTSEC,ZZ 
COMMON/SOLID/NBDIES,NWALL(2),WALL(215,l), 
>ZCR(215,1),HUB(2),Z0(2),THETA(215,1),NINC(2), 

>INC(15,2),XMAX(2) 

COMMON/SYSTEM/NDIM,A(215,215),X(215),IPVT(215) 

DEFINE SOLID. 

CALL SOLID(NBDIES, NWALL, WALL, CHARD) 

COMPUTE ARCLENGTH. 

ARCL=0 
NDIM=0 

DO © L=l, NBDIES 
NDIM=NDIM+NWALL(L) 

DO © K=l, NWALL (L) 

0 ARCL=ARCL+CABS(WALL(K,L)-WALL(l+MOD(K,NWALL(L)),L)) 

C 

C COMPUTE CORE RADIUS. 

R0=ARCL/NDIM/2. 
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SIGMA=R0/3 
WRITE(6, 104) SIGMA 

104 FORMAT(/," THE CORE RADIUS WAS CHOSEN AS: '\F8.4,/) 
SIGMA2=SIGMA**2 

DEFINE CREATION POINTS. 

D06 L=l,NBDIES 
DOT K=1,NWALL(L) 

ZZ=WALL(l+MOD(K+NWALL(L)-2,NWALL(L)),L) 

> -WALL(l+MOD(K,NWALL(L)),L) 

7 ZCR(K,L)=WALL(K,L)+CMPLX(0.,R0/CABS(ZZ))*ZZ 
C CHECK THAT WALL AND CREATION POINTS ARE NOT 
C CROSSED DUE TO TOO SHARP A CONCAVE KINK OR AN 
C ERROR IN DEFINING THE BODY. 

XMAX(L)=REAL(ZCR(1,L)) 

D06 K=1,NWALL(L) 

XMAX(L)=AMAXl(XMAX(L),REAL(ZCR(K,L))) 

KP=l+MOD(K,NWALL(L)) 

IF(REAL((ZCR(KP,L)-ZCR(KD))* 

> CONJG(WALL(KP,L)-WALL(K,L)))GT.O.)GO TO 6 
WRITE (6,102)L,K 

102 FORMAT ( 

>” ON BODY NUMBER ",13," YOU HAVE TOO SHARP A", 

>” KINK NEAR POINT " 44) 

STOP 

6 CONTINUE 
C 

C COMPUTE MATRIX. 

C EXCEPT FOR THE LAST ONE ON EACH BODY, EACH LINE 
C WILL REPRESENT: PSI(WALL(NWALL))-PSI(WALL(I)) 

C WHERE PSI CORRESPONDS TO THE NEW VORTICES TO BE 
C CREATED AND WILL HAVE TO CANCEL THE PSI DUE TO 
C THE FREESTREAM + OLD VORTICES. 

C THE LAST LINE IS ALL 0, EXCEPT ON COLUMNS 
C BELONGING TO THE SAME BODY, THEN IT IS 1. 

C IT WILL CONTROL THE TOTAL STRENGTH OF ALL THE NEW 
C VORTICES EMANATING FROM THAT BODY, AS WELL AS 
C PREVENTING THE MATRIX A FROM BEING SINGULAR. 
PI=4*ATAN(1.) 

10=0 

D02 Ll = l,NBDIES 
J0=0 

D04 L2»1,NBDIES 
KRON=0 
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IF(Ll.EQ.L2)KRON=l 
D03 J=1,NWALL(L2) 

A(IO+NWALL(L1),JO+J)=KRON 
D03 I=1,NWALL(L1)-1 

3 A(IO+I,JO+J)=-.25/PI*ALOG( 

> (SIGMA2+CABS(WALL(I+1,L1)-ZCR(J,L2))**2) 

> /(SIGMA2+CABS(WALL(I,L1)-ZCR(J,L2))**2)) 

4 JO=JO+NWALL(L2) 

2 I0=I0+NWALL(L1) 

C 

C GAUSS ELIMINATION. 

C SGECO IS THE LINPACK ROUTINE. 

C X IS USED AS A DUMMY HERE. 

CALL SGEC 0( A,2 1 5,NDIM,IPVT ,COND ,X) 

COND=l/COND 
WRITE (6,103)COND 

103 FORMAT(/," CONDITION NUMBER " ,E8.2) 

C 

C MISCELLANEOUS; BISECTION DEVICE. 

DOl L=1,NBDIES 

C FIND A HUB FOR THE RAYS OF THE BISECTION DEVICE. 
HUB(L)— 0 
NPTS=0 
FERST=0 

DO 20 I=1,NWALL(L) 

EP=l+MOD(I,NWALL(L)) 
IM=l-f-MOD(I+NWALL(L)-2,NWALL(L)I 
IF(1 ,+AIMAG((ZCR(IP,L)-ZCR(I,L))* 

> CONJG(ZCR(IMT-)-ZCR(I ) L))).EQ.l.)GO TO 20 
NEXT(I)=FERST 

FIRST=I 

20 CONTINUE 
C 

I*=FIRST 

21 IF(I.EQ.O)GO TO 22 
IP=l-fMOD(I > NWALL(L)) 

K=NEXT(I) 

24 IF(KEQ.O)GO TO 23 

KP=l+MOD(K,NWALL(L)) 

AA=AIMAG( (ZCR(IP,L)-ZCR(I,L)) * 

> CONJG (ZCR(KP,L)-ZCR(K,L)) ) 

IF(AA EQ.O.)GO TO 10 

INTSEC=ZCR(I,L)+AIMAG((ZCR(K,L)-ZCR(I,L))‘ 

> CONJG(ZCR(KP,L)-ZCR(K,L)))* 
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> (ZCR(IP,L)-ZCR(I,L))/AA 
DO 11 J==1,NWALL(L) 

11 IF(1.+AJMAG ( (INTSEC-ZCR(J,L)) * 

> CONJG(ZCR(l+MOD(J,NWALL(L)),L)-ZCR(J,L)) 

> ) .LT. l.)GO TO 10 
HUB(L)=HUB(L)+INTSEC 
NPTS=NPTS-j- 1 

10 K=NEXT(K) 

GO TO 24 
23 I s= NEXT (I) 

GO TO 21 

22 IF(NPTS.GE.l)GO TO 12 
WRITE(6,100)L 

100 FORMAT { 

>" THERE IS A PROBLEM WITH BODY NUMBER " ,14,/, 

>" ITS SHAPE IS PROBABLY TOO COMPLEX AND THE " 

>, PROGRAM WAS UNABLE TO DEFINE A HUB. " , 

>" OR ELSE YOUR POINTS ARE NOT COUNTERCLOCKWISE”) 
STOP 

12 HUB(L)=HUB(L)/NPTS 

WRITE (6,101)L,REAL(HUB(L))AIMAG(HUB(L)) 

101 FORMAT (/," THE HUB FOR BODY * ,14,” IS AT: " ,2F8.3) 

C NOW COMPUTE THE POLAR ANGLE OF ALL THE POINTS. 
ZO(L)=CONJG(ZCR(NWALL(L),L)-HUB(L)) 

D08 1=1 ,NWALL(L) 

8 THETA(I,L)= 

> AIMAG(CLOG(-Z0(L)*(ZCR(I,L)-HUB(L)))) 

c 

C NOW COMPUTE THE INCREMENTS FOR THE BISECTION. 
INC(l,L)=NWALL(L)/2 
K=NWALL(L)-INC(1 JL) 

DOS 11=2,15 
EF(K.EQ.l)GO TO 1 
INC(II,L)=MAX 0 ( 1 ,INC(n- 1 ,L)/2) 

5 K=K*INC(II,L) 

1 NINC(L)=n-l 
C 

RETURN 

END 

SUBROUTINE SOLID(NBDIES,NWALL, WALL, CHARD) 

C 

C ALLOWS THE USER TO INPUT THE SHAPE OF THE SOLID. 

C IT CAN BE 1 BODY (NBDIES=1), OR SEVERAL. 

C NWALL(L) NUMBER OF POINTS ON SOLID L (L=1,NBDIES). 
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C WaLL(K,L) complex positions of these POINTS 
C (K=1,NWALL(L)). 

C THEY MUST RUN COUNTER-CLOCKWISE!! 

C CHARD IS THE CHARACTERISTIC DIMENSION 
C 

COMPLEX WALL(215,1) 

INTEGER NW\LL(2) 

C 

NBDIES=’ 

C FIRST BODY. 

NWALL(1)=200 
DO 1 K=l,100 
Y=-1+.02*K 

WALL(K,1)=CMPLX(.2*(1-Y*Y),Y) 

1 WALL(K+100,l)=CMPLX(0.,-Y) 

CHARD =2 
C 

WRITE(6,100) 

100 FORMAT (/," THE SOLID IS A CAMBERED FLAT PLATE") 
C 

RETURN 
END 

SUBROUTINE INIT(ISTART,NSTART,T,VO,ABSUIN, 

> CHARD, DO) 

INITIALIZE TIME-DEPENDENT VARIABLES. 

COMPLEX Z,V,VM 

COMMON/VORTEX/NVORT,Z(2000),V(2000),VM(2G00) 

> ,GAMMA(2000) 

IF(ISTART.EQ.O)GO TO 8 

CASE OF A START FROM POTENTIAL FLOW. 

START AT STEP 1 WITH T=0 AND NO VORTICES. 

NS TART =1 
T=0 

NVORT =0 

C GIVE PHONY VALUES TO THE VORTEX POSITIONS AND 
C CIRCULATIONS. 

D07 1=1,2000 
Z(I)=10. 

VM(I)=0. 

7 GAMMA(I)=0. 

C A TENTATIVE VALUE FOR V0, WHICH WILL BE 
C ADJUSTED LATER. 
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V0 = 1.E-5*ABSUIN*(CHARD/D0)**3 
RETURN 

CASE OF A START FROM A PREVIOUS RUN. 

READ (9) NSTART,T ( NVORT,Z,GAMMA,VM ( VO 

RETURN 
END 

SUBROUTINE MERGEIDO.VO.NDES.N) 

THIS ROUTINE MERGES PAIRS OF VORTICES INTO ONE 
WHENEVER THE PENALTY FOR DOING SO 
IS UNDER A TOLERANCE (VO). 

ALSO ADJUSTS VO TO ACHIEVE DESIRED 
NUMBER OF VORTICES. 

REAL B(2000),MERTST(2000) 

COMPLEX HUB.Z.ZI/V.ZCR.ZO.WALL.VM 
COMMON/ VORTEX/NVOR7 ,Z(2000),V(2000),VM(2000) 

> ,GAMMA(2000) 

COMMON/SOLED/NBDIES,NWALL(2),WALL(215,l), 

> ZCR(215,I)TIUB(2),Z0(2) ) THETA(215,1),NINC(2) 

> 4NC(15.2),XMAX(2) 

RETURN IF THERE ARE NO VORTICES YET. 
IF(NVORT.LE.5)RETURN 

FEEDBACK NUMBER OF VORTICES TO THE TOLERANCE. 
V0=V0*EXP(AMAX1(-.1,.00<.*(NV< >RT-NDES))) 
SQV0=SQRT(V0) 

PREPARE VORTICES FOR MERGING; 

J FIND THEIR DISTANCE TO THE WALL. 

DO 7 I~l,NVORT 
D*l.E+30 
DOl L=I,NBDIES 

C FIND THE PROJECTION ON EACH WALL BY BISECTION 
C OF THE POLAR ANGLE. 

TTA*.AIMAG(CLOG(-ZO(L)*(Z(I)-HUB(L)))) 

K=NWALL(L) 

DOS n=l,NINC(L) 

5 IF(T1A.LT.THETA(K-INC(II,L),L)) 

> K»K-INC(1I,L) 

KM=l-f-MOD(K-+-NW'ALL(L)*2,NW.ALL(L)) 

C D IS THE DISTANCE TO THE WALL. B(I) IS STORED 
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C TO REDUCE THE WORK IN THE INNER LOOP 
C (LABEL 4, SEE FURTHER DOv/N) TO A MINIMUM. 

1 D=AMINl(D,AIMAG(CONJG(Z(I)-ZCR(K,L))* 

> (ZCR(K,L)-ZCR(KM,L)))/CABS(ZCR(K,L)-ZCRpvM,L))) 
B(I)=ABS(GAMMA(I))*(ABS(D+D0))**(-1.5)/SQV0 

7 CONTINUE 
C 

C TAKE THE VORTICES ONE BY ONE, STARTING WITH 
C THE LAST ONE. 

NVORTO=NVORT 
DO 3 I=NVORTO,2,-1 

TEST ALL THE OTHER VORTICES AGAINST IT. 

DO 4 J=1,I-1 

MERTST(J)=(HEAL(Z(J)-Z(I))**2+AIMAG(Z(J)-Z(I))**2)* 

> B(J)*B(I)-ABS(GAMMA(J)+GAMMA(I)) 

IS THE MERGING TOLERANCE SATISFIED? 

(ISMIN IS A CRAY FUNCTION USED HERE TO FIND THE 
MINIMUM OF "MERT5T” 

FROM 1 TO IM WITH INCREMENT 1.) 

J=ISMIN(I-1 ,MERTST, 1) 

IF NOT SO GO TO NEXT INDEX I. 

IF(MERTST(J).GT.O.)GO TO 3 
IF SO, PROCEED WITH THE MERGING: 

PUT THE NEW VORTEX IN J’S PLACE. 

Z(J)- (Z(I)*GAMMA(I)+ Z(J)*GAMMA(J))/ 

> (GAMMA(I)+GAMMA(J)) 

VM( J) = ( VM (I) * GAMMA(I ) + VM( J) * G AMM A( J) ) / 

> (GAMMA(I) + GAMMA( J)) 
GAMMA(J)=GAMMA(I)+GAMMA(J) 
B(J)=B(I)*ABS(GAMMA(J)/GAMMA(I)) 

PUT LAST VORTEX (INDEX: NVORT) IN Ith PLACE. 
Z(I)=Z(NVORT) 

VM(I)= VM(NVORT) 

GAMMA(I)=GAMMA(NVORT) 

B(I)=B(NVORT) 

NVORT =NVORT-l 

TAKE NEXT VORTEX ON THE LIST. 

3 CONTINUE 
C 

END 


102 


o o o o o 


OP Pr» 

■ * * w 


*> V'-.’ii.i’iY 


SUBROUTINE BCBODY(FORCE, MOM, GAMMAO,N,T,UINF, 

> CHARD, NOLD, CP, DELT,SIGMA2) 

COMPLEX FORCE(2),UINF 

REAL MOM(2),DPDS(215,l),CP(215,l) 

C 

C THE BODY ABSORBS VORTICES AND EMITS NEW VORTICES 
C TO ACCOUNT FOR IT, PLUS SOME NEW VORTICITY WHICH 
C WILL ALLOW THE VELOCITY FIELD 
C TO SATISFY THE BOUNDARY CONDITION U=V=0. 

C 

C DETECT AND ABSORB VORTICES THAT CRASHED INTO WALL. 
C START COMPUTING PRESSURE AND FORCE. 

CALL ABSORB(DPDS, FORCE, MOM, GAMMA0,N) 

C 

C EMIT NEW VORTICES TO SATISFY BOUNDARY CONDITION, 

C AND FINISH COMPUTING PRESSURE, FORCE, ETC. 

CALL EMIT(N,T,UINF,CHARD,NOLD, FORCE, MOM, DPDS, 

> CP,DEIT,SIGMA2) 

C 

RETURN 
END 

SUBROUTINE ABSORB(DPDS, FORCE, MOM.GAMMAO.N) 

REAL DPDS(215,l),MOM(2) 

COMPLEX FORCE(2),Z,V,VM, WALL, ZCR, HUB, ZO 
COMMON/VORTEX/ NVORT,Z(2U00),V(2000),VM(2000) 

> ,GAMMA(2000) 

COMMON/SOLID/NBDIES,NWALL(2),WALL(21b,l), 
>ZCR(215,1),HUB(2),Z0(2),THETA(215,1),NINC(2) 
>,INC(15,2),XMAX(2) 

COMMON/SYSTEM/NDIM,A(215,215),X(215),IPVT(215) 

KILL THE VORTICES THAT ARE TOO CLOSE TO A WAI L. 

(LOST VORTICITY WILL BE REINTRODUCED IMMEDIATELY.) 

TAKE THE SEPARATE BODIES ONE BY ONE. 

10=0 

D09L=1,NBDIES 
C GET READY TO COMPUTE THE FORCE, MOMENT, AND \V.\LL 
C PRESSURE ON BODY "L". 

FORCE(L)=0 

MOM(L)=0 

DO 3 K=1,NWALL(L) 

3 DPDS(K,L)=0 
I0=I0+NWALL(L) 
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C X(IO) IN THE SECOND MEMBER WILL BE THE CIRCULATION 
C DEFECT TO BE FILLED UP BY THE NEW VORTICES 
C EMANATING FROM THIS BODY. 

X(I0)=0 

C IF THIS IS THE FIRST STEP, INTRODUCE THE DESIRED 
C CIRCULATION AROUND THE FIF.ST BODY. 
IF(N.EQ.1.AND.L.EQ.1)X(I0)=GAMMA0 

LOOK AT THE VORTICES ONE BY ONE. 

1=1 

IF(I.GT.NVORT)GO TO 9 
IF(REAL(Z(I)).GTXMAX(L))GO TO 5 
OVER WHICH WALL SEGMENT IS THE VORTEX? 

FIND IT BY BISECTION OF THE POLAR ANGLE. 
TTA=AIMAG(CLOG(-ZO(L)*(Z(I)-HUB(L)))) 

K=NWALL(L) 

DO 7 n=l,NINC(L) 

IF(TTA.LT.THETA(K-INC(n,L),L)) 

> K=K-INC(II,L) 

MAKE SURE OF WHERE THE PROJECTION OF THE VORTEX 
IS ON THE WALL. 

KM=l+M0D(K-2+NWALL(L),NWALL(L)) 
IF(REAL((Z(I)-ZCR(KM,L))*CONJG(ZCR(K,L)- 

> ZCR(H M0D(K+NWALL(L)-3,NWALL(L)),L))).GE 0 K<0 TO 1 
K=KM 
GO TO 2 

KP=1+M0D(K,NWALL(L)) 
IF(REAL((Z(I)-ZCR(K,L))*CONJG(ZCR(KP ■ , 

> •ZCR(KM,L)))EE.O.)GO TO 4 
KM=K 
K=KP 
GO TO 1 

IS THE VORTEX INSIDE THE S^lID? 

KM=l+M0D(K-2+NV .L,(L),NWALL(L)) 
D=AIMAG(C0NJG(7' oCR(K,L))*(ZCR(K,L)-ZCR(KM,L))) 

IF IT IS NOT, LEA\~ T ALONE. 

IF(D.GT.O )G f '35 
C IF IT IS, K3LI » \ iFIRST RECORD THE LOSS OF 
C CERCULATI AND LINEAR AND ANGUL AR MOMENTUM) 

X(IO' - aJO)+GAMMA(I) 

FO E(L)=FORCE(L)-GAMMA(I)*Z(I) 
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MOM(L)=MOM(L)-GAMMA(l)*(REAL(Z(I)-HUB(L)) 

> **2-f-AIMAG(Z(I)-HUB(L))**2) 

C ALSO RECORD VORTICITY ABSORPTION FOR THE PRESSURE 
C GRADIENT. LINEAR DEPOSITION. 

D=REAL((Z(I)-ZCR(KM,L))* 

> CONJG(ZCR(K,L)-ZCR(KM,L))) 

DM=REAL((ZCR(K,L)-Z(l))* 

> CONJG(ZCR(K,L)-ZCR(KM,L))) 
DPDS(K,L)=DPDS(K,L)+GAMMA(I)*D/(D+DM) 
DPDS(KM,L)=DPDS(KM,L)+GAMMA(I)*DM/(D+DM) 

NOW PUT LAST VORTEX IN THE Itb PLACE IN THE ARRAY. 
Z(I)=Z(NVORT) 

GAMMA(I)=GAMMA(NVORT) 

VM(I)=VM(NVORT) 

NVORT=NVORT-l 
GO TO 6 

GO TO NEXT VORTEX. 

5 1=1+1 

GO TO 6 
0 CONTINUE 
C 

RETURN 

END 

SUBROUTINE EM1T{N,T,U1NF,CHARD,N0LD, FORCE, 

> MOM,DPDS,CP,DELT,SIGMA2) 

COMPLEX Z.V/VM, WALL, ZCR, HUB, ZO,FORCE(2),UINF 
REAL MOM(2),DPDS(215,l),CP(215.1),PSI(215,l) 

> T*S(2000) 

CCMMON/VORTEX/NVORT,Z(2000),V(2000),VM(2()00) 

> ,GAMMA(2000) 

COMMON/SOLH)/NBDIES,NWALL(2),WALL(215,l), 

> ZCR(215,1),HUB(2),Z0(2),THETA(215,1),NINC(2) 

> 4NC(15,2),XMAX(2) 

COMMON/SYSTEM/ND IM , A(2 15, 215) ,X(2 15), IPVT(2 15) 

C 

PI=4*ATAN(1 .) 

10=0 

DOS L RNBD1ES 

C STREAM FUNCTION AT THE WALL; CONTRIBUTION OF... 

DO lfl K«1,NWALL(L) 

C ...THE FREESTREAM... 

C ...AND OF THE OLD VORTICES. 

D04 I-l.NVORT 
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4 PS(I)=GAMMA(I)* 

> ALOG(REAL(Z(I)-WALL(K,L))**2+ 

> AIMAG(Z(I)-WALL(K,L))**2-)-SIGMA2) 

16 PSI(K,L)=AIMAG(WALL(K,L)*CONJG(UINF)) 

> -SSUM(NVORT,PS,l)*. 25/PI 
C 

C THE NEW VORTICES MUST CANCEL THE STREAM FUNCTION 
C AT THE WALL. 

C THIS GIVES A LINEAR SYSTEM FOR THEIR CIRCULATIONS. 

C COMPUTE ITS SECOND MEMBER. 

DOl8 K=1 ,NWALL(L)- 1 
18 X(I0-H<)=PSI(K,L)-FSI(K+1,L) 

5 I0=I0+NWALL(L) 

C 

C SOLVE SYSTEM. 

SGESL IS THE LINPACK ROUTINE. 

X IS THE RIGHT HAND SIDE, AND 
THE SOLUTION IS WRITTEN OVER IT. 

CALL SGESL(A,215,NBODY,IPVT ,X,0) 

CREATE NEW VORTICES. 

C REMEMBER WHICH VORTICES ARE OLD ENOUGH TO USE 
C ADAMS-BASHFORTH. 

NOLD=NVORT 

10=0 

DO 17 L=1,NBDIES 
DO 3 K=1,NWALL(L) 

C PUT THE NEW VORTEX AT THE END OF THE ARRAY. 

NVORT=NVORT+l 
Z(NV ORT ) =ZCR(K,L) 

GAMMA(NVORT) =X(I0 +K) 

C RECORD THE GAIN OF LINEAR AND ANGULAR MOMENTUM. 
FORCE(L)=FORCE(L)+GAMMA(NVORT)*Z(NVORT) 
MOM(L)=MOM(L)+GAMMA(NVORT)*(REAL(Z(NVORT) 

> -HUB(L))**2+A1MAG(Z(NV0RT)-HUB(L))*‘2) 

C ALSO RECORD VORTICITY CREATION FOR THE PRESSURE 
C GRADIENT. 

3 DPDS(K,L)=DPDS(K,L)-X(IO+K) 

C 

C FILTER PRESSURE GRADIENT AND INTEGRATE IT TO GET 
C PRESSURE 

CP(1,L)=0 | 

D014 K=2,NWALL(L) \ 

i 
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CP(K,L)=CP(K-l,L)+(3*(DPDS(K,L)+DPDS(Iv-l,L)j 

> +DPDS(l+MOD(K+NWALL(L)-3,NWALL(L>),L)+ 

> DPDS(l+MOD(K,NWALL(L)),L))/ 

> (4*DELT*CABS(UINF)**2) 

C 1RMALIZE PRESSURE. 

CPMAX=CP(ISMAX(NWALL(L),CP(l,L),l),L) 

D022 K=1,NWALL(L) 

22 CP (K,L) = CP(K,L)- CPMAX 

C 

C FINISH COMPUTING FORCE AND MOMENT AND 
C NON-DIMENSIONALIZE THEM. 

FORCE(L)=FORCE(L)* 

> CMPLX(0.,2./(DELT*CHARD*CABS(UINF)**2)) 
MOM(L)=MOM(L)*2/(DELT*(CHARD*CABS(UINF))**2) 

17 I0=I0+NW 
C 

C PRINT INSTANTANEOUS DATA. 

T=T+DELT 

IF(MOD(N,5)£Q.O) 

> WRITE (6,105)N,T,NVORT,REAL(FORCE(l)), 

> AIMAG(FORCE(l)),MOM(l) 

105 FORMAT (/,6H STEP ,I4,6H TIME ,F8.4,7H NVORT , 

>14, 3H CD,F7.4,3H CL, 

>F7.4,4H MOM,F7.4) 
IF(MOD(N,5).EQ.OAND.NBDIES.GE.2) 

> WRITE(6 , 106)REAL(FORCE( 2)),AIMAG(FORCE(2) ) 

> ^!OM(2) 

106 FORMAT (” CD, CL, AND MOM ON SECOND BODY:” 
>,2F10.4,F1 1.4) 

RETURN 

END 

SUBROUTINE MOVE) UINF,SIGMAi ,N OLD ,DELT ) 

C 

C MOTION OF THE VORTICES. 

C OLD VORTICES USE ADAMS-BASHFORTH-2, 

C NEW ONES USE EULER EXPLICIT. 

C 

COMPLEX Z,V,VM,UINF 

COMMON/VORTEX/NVORT,Z(2000),V(2000),VM(2000) 

> ,GAMMA(2000) 

C 

C COMPUTE VELOCITIES OF THE VORTICES. 

CALL VELOC T (UINF.SIGMA2) 

C 

C MOVE VORTICES. 

C ADAMS-BASHFORTH-2 FOR THE OLD VORTICES. 
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D013 1=1, MOLD 

Z(I)=Z(I)+DELT*(1.5*V(I)-.5*VM(I)) 

13 VM(I)=V(I) 

C EULER EXPLICIT FOR THE NEW VORTICES. 

DO 2 I=NOLD+l,NVORT 
Z(I)=Z(I) -j-DELT *V (I) 

2 VM(I)=V(I) 

C 

RETURN 
END 

SUBROUTINE VELOCT(UINF,RC2) 

BIOT SAVART INTERACTION OF VORTICES, POSITIONS Z(I), 
CIRCULATION GAMMA(I). VELOCITY AT INFINITY =UINF . 

PC IS THE CHARACTERISTIC RADIUS IN THE CUT-OFF: 
U(R) IS= (GAMMA/2PI)*R/(R**2+RC**2) 

REAL VX(2000),VY(2000) 

COMPLEX Z,V,DELZ ,VM,UINF 

COMMON/VORTEX/ NVORT,Z( 2000), V(2000),VM(2000) 

> ,GAMMA(2000) 

FREESTREAM VELOCITY. 

PI=4*ATAN(1.) 

DO 3 l=l,NVORT 
V(I)=CMPLX(0.,-2*PI)*UINF 

COMPUTE INTERACTIONS. 

LOOP ON FIRST VORTEX. 

DOl I=2,NVORT 
LOOP ON SECOND VORTEX. 

D04 J= 1,1-1 
DELZ=Z(I)-Z(J) 

DELZ=DELZ/(RC2+REAL(DELZ)**2+AIMAG(DELZ)**2) 
VX( J)=REAL (DEL Z)* GAMM A( J) 
VY(J)=AIMAG(DELZ)*GAMMA(J) 
V(J)=V(J)-GAMMA(I)*DELZ 

(THE CRAY FUNCTION SSUM IS USED TO SUM UP A VECTOR 
LIKE VX, FROM 1 TO 1-1 WITH INCREMENT 1.) 

(NOTE THAT THE CSUM FUNCTION WOULD HAVE BEEN THE 
LOGICAL CHOICE HERE; BUT IT SEEMS TO HAVE A BUG...) 
V(I)=V(I)+CMPLX(SSUM(I-1,VX,1),SSUM(I-1,VY,1)) 

MULTIPLY BY I/2PI. 

D02 I=*l,NVORT 
2 V(I)=aV(I)*CMPLX(0. ,.5/PI) 


RETURN 

END 
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Fig. 5. Division of the flow field into an inviscid and a viscous domain 
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Fig. 6b. Core Velocity. Point Vortex, Core 1, 
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Fig. 6c. Core Stream Function. Point Vortex, 

Core 1, •••Core 2 
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Fig. 9. Schematic of cell-to— cell interactions 



Fig. 10a. Same sign. 



Fig. 10. Schematic of two cases of merging 
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Fig. 13. Displacement of tut* outer flow by the Vortex Sheet. 



Fig t4. Flow Chart of lYogrnm KPDH 
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Fig. l6a. Simulation with 0"=O. 005 
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Fig. 17b. Simulation with 5 =.04 
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Fig. 18a. Simulation with £=0 



Fig. 18b. Simulation with £ =0.8 10 
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Fig. 24. Stills of the Dynamic Stall Simulation 
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Fig. 26c. Out-of-phase pressure 
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Fig. 27. Results of the tilt-rotor study. 
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Fig. 31. Shedding frequency of circular cylinder. 
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Angle 9 from front of cylinder 

Fig. 33a. Pressure Distribution. Exp. , Re=10 6 [62] 

Comp. , Re=10 4 
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Fig. 33b. Pressure distribution. Exp., Re=10 8 [62] 
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Fig. 33c. Pressure Distribution. Exp., Re~2.6 10* [63] 

Comp., Re=3.16 10 8 
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Angle 9 from front of cylinder 

Fig. 33e. Pressure Distribution. Exp., Re=8.5 10* [63] 

Low resolution, Re=10* • • -High resolution, Re=10* 


Angle 9 from front of cylinder 

Fig. 33f. Pressure Distribution. Exp., Re=3.6 10* [63] 
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Fig. 33g. Pressure Distribution. Exp., Re=8.5 10* [64] 

Comp., Re=10 7 
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Fig. 34b. Friction Distribution. Exp., Re=10*[63j 
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Fig. 34c. Friction Distribution. Exp., Re=2.6 10 r [63] 

- - Comp., Re=3 16 10 6 
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Fig. 34e. Friction Distribution. Exp., Re=8.5 10®[63] 

Low resolution, Re=10® High resolution, Re=10® 



j Angle 9 from front of cylinder 

1 Fig. 34f. Friction Distribution. Exp., Re=3.6 lO'fCS] 

Comp., Re=3.16 10® 


145 


i 



Fig. 34g. Friction Distribution. Comp., Re-10 r 
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Fig. 35b. Simulation at Re=10 6 



Fig. 38. Simula:. 
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Fig. 38. Velocity on centerline Exp., Re=1.4 10 s [62] 

Computations, Re=10 9 


149 







Amplitude Amplitude Amplitude 

0.0 0.1 0.2 0.3 0.4 0.6 0.8 0.0 0.1 0.2 0.3 0.4 0.6 0.6 0.0 0.1 0.2 0.3 0.4 


original PA ®f 5 
OF POOR QUALITY 



Fig. 41a. Re=10 9 



Fig. 41b. Re=3.16 10 s 



Fig. 41c. Re=3.16 10® 
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Fig. 44. Detail near wall in programs KPD 1 and X PD 1. 



Fig. 45. Flow Chart of Program KPD1 
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Fig. 46. Flow Chart of Program KPD2 
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